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ABSTRACT 


Analytic  descriptions  of  the  transient  thermal  field 
are  presented  for  three  configurations  of  electrodes  carry¬ 
ing  direct  current;  two  concentric  spheres,  a  single  sphere 
in  an  infinite  medium,  and  two  infinitely  long  co-axial  cy¬ 
linders.  The  medium  characteristics  are  taken  to  be  indepen¬ 
dent  of  temperature  and  electric  field.  Homogeneity  and  iso¬ 
tropy  are  also  assumed. 

The  electric  boundary  conditions  used  are  constant 
potential  at  the  inner  shell,  zero  potential  at  the  outer 
(or  infinity) .  The  thermal  solution  is  derived  on  the  basis 
of  no  heat  flow  across  the  inner  shell  surface  and  constant 
(but  arbitrary)  temperature  at  the  outer. 

Sample  computer  programs  are  given  for  the  first 
two  cases. 
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INTRODUCTION 

With  the  development  of  electrical  demand  has  come 
the  economic  feasibility  of  high  voltage  direct  current  trans¬ 
mission  over  relatively  large  distances.  The  reliability  and 
economy  of  these  circuits  may  be  greatly  enhanced  by  the  per¬ 
manent  or  emergency  use  of  the  earth  as  a  return  conductor, 
enabling  the  D.C.  link  to  operate  at  half  capacity  in  the 
event  of  the  failure  of  one  D.C.  pole.  The  use  of  the  earth 
as  a  return  circuit,  however,  necessitates  making  stable 
ground  connections  at  the  terminals.  It  is  one  of  the  prob¬ 
lems  associated  with  these  connections  with  which  this  dis¬ 
sertation  is  concerned. 

Examination  of  reports  concerning  projects  in  high 
voltage  D.C.  transmission  in  the  western  world,  viz.  the  Got¬ 
land  link  in  Sweden,  the  Bonneville  Power  Administration  link 
in  the  western  United  States,  and  the  New  Zealand  tie,  as  well 
as  papers  on  electro-osmosis,  revealed  essentially  two  areas 
open  to  immediate  further  development.  The  first  of  these 
areas  concerns  the  problems  of  relating  thermal  and  electri¬ 
cal  properties  of  the  soil  to  electric  field,  temperature, 
and  soil  properties  such  as  grain  size,  density,  material  com¬ 
position  and  moisture  content.  The  second,  perhaps  of  more 
urgent  practical  interest,  is  that  of  transient  thermal  be¬ 
havior  of  ground  electrodes. 
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In  the  specific  area  of  thermal  behavior,  theory  has 

3 

been  developed  for  two  restrictive  conditions.  Rudenberg 

gives  long-term  (steady  state)  solutions  of  temperature  fields 

surrounding  various  electrodes  and  predicts  a  "time  constant" 

on  the  basis  of  steady  rise  of  temperature  at  the  initial  rate. 

4  5  11  12  13 

Other  authors  '  '  '  '  either  quote  these  solutions  or  de¬ 

velop  the  same  conclusions  independently. 

In  this  thesis  the  restrictions  on  time  are  completely 
removed.  Three  configurations  of  electrodes  are  considered. 

The  first  involves  two  concentric  spheres  with  heat  generation 
between  them,  the  second  is  a  single  sphere  embedded  in  an  in¬ 
finite  medium,  and  the  third  consists  of  two  co-axial,  infin¬ 
itely  long  cylinders.  The  electrical  and  thermal  performance 
of  these  three  arrangements  is  analysed  and  formal  solutions 
are  presented.  Computer  programs  and  particular  numerical 
solutions  are  included  for  the  first  two  electrode  arrange¬ 
ments  . 


In  deriving  the  thermal  field  solutions  some  mathe¬ 
matical  complexity  has  been  avoided  by  assuming  that  the 
electrical  and  thermal  parameters  remain  constant,  independent 
of  temperature  or  electric  potential.  The  medium  is  further 
assumed  to  be  homogeneous  and  isotropic.  In  actual  fact  soil 
has  a  negative  temperature  coefficient  of  electrical  resist¬ 
ivity  which  would  make  these  results  slightly  pessimistic. 

On  the  other  hand,  the  assumptions  of  homogeneity  and  isotropy 
are  certainly  open  to  question,  and  could  be  important  consider- 
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ations  if  they  tended  to  cause  low  predictions  of  temperature. 

For  the  arrangements  with  two  concentric  spheres  or 
cylinders  the  boundary  conditions  used  are  as  follows: 

1)  the  inner  shell  is  maintained  at  constant  D.C. 
potential 

2)  the  outer  shell  is  held  at  reference  or  zero 
potential 

3)  there  is  no  heat  flow  across  the  inner  shell  sur¬ 
face,  that  is  the  thermal  gradient  at  the  inner 
shell  is  always  zero 

4)  the  outer  shell  is  held  at  its  initial  temperature 
for  all  time. 

Conditions  on  the  single  sphere  problem  are  essenti¬ 
ally  the  same  if  infinity  is  considered  as  the  outer  shell. 

Condition  1  above  closely  approximates  a  practical 
situation  providing  the  assumptions  regarding  the  thermal 
and  electrical  parameters  of  the  medium  are  realistic.  In 
this  case  a  constant-current  source  (practical  electrode) 
and  a  constant  potential  source  are  identical  as  far  as  analy¬ 
sis  is  concerned.  The  equivalence  of  the  duals  is  established 
more  rigorously  in  section  A-2,  and  the  entire  development 
may  be  carried  out  on  the  assumption  of  constant  current  if 
so  desired. 

The  electrodes  themselves  will  normally  be  constructed 
of  materials  which  are  relatively  good  electrical  conductors. 
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For  this  reason  very  little  heat  will  be  generated  within  the 
electrode  itself  in  comparison  with  that  produced  between  the 
two  shells.  Assumption  3  is  equivalent  to  saying  that  the 
amount  of  heat  generated  within  the  inner  electrode  is  just 
sufficient  to  match  its  temperature  with  that  of  the  immediate 
surroundings.  In  actual  practice  this  may  or  may  not  be 
strictly  true  but  should  be  realistic  enough  to  render  the 
solutions  useful.  If  heat  actually  flows  into  the  electrode 
the  solutions  presented  will  be  slightly  pessimistic  in  that 
they  will  predict  a  temperature  higher  than  that  actually  ob¬ 
served  . 

The  concentric  shell  problems  were  initially  attempted 
as  an  approach  to  the  infinite  medium  case.  Assumptions  2  and 
4  were  therefore  based  on  practical  considerations  from  the 
outset;  electric  potential  is  invariably  measured  with  res¬ 
pect  to  a  reference  far  removed  from  the  electrode,  and  at 
large  distances  from  the  electrode  heat  generation  in  the 
ground  will  be  negligible. 

The  medium  is  initially  assumed  to  be  at  a  uniform 
temperature  equal  to  ambient.  Solutions  for  other  initial 
conditions  are  possible,  at  least  in  principle,  as  discussed 
in  section  5.2.  Transient  electric  behavior  is  assumed  suf¬ 
ficiently  rapid  as  to  be  instantaneous  in  comparison  with 
thermal  time  dependence.  The  initial  electric  condition  is 
therefore  a  step  rise  in  potential  (or  current)  from  zero  to 


its  final  value. 
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The  solutions  given  here  assume  complete  radial  sym¬ 
metry.  However,  they  are  reducible  to  the  more  realistic 
case  of  a  half  shell  buried  so  as  to  have  its  flat  surface 
flush  with  the  ground  by  merely  doubling  the  electrode  cur¬ 
rent  used  in  the  computations. 

In  order  to  avoid  needless  repetition  the  symbols 
used  in  the  following  analyses  and  in  the  appendices  are  ' 
listed  on  the  following  pages  and  are  not  re-defined  when 
they  appear  in  the  text. 

References  are  indicated  by  a  superscript  which  refers 
to  the  numbering  of  the  bibliography  entries.  Parentheses 
around  an  equation  number  indicate  repetition  of  an  equation 
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Electrode  Arrangement  for 
Concentric  Sphere  Problem 


1.1 


CHAPTER  1:  CONCENTRIC  SPHERES 

1 . 1  Introduction 

This  chapter  deals  with  the  analytic  solution  for 
the  electric  and  thermal  fields  in  a  region  between  two  con¬ 
centric  spherical  shells  (see  fig.  1.1).  As  indicated  in  the 
"Introduction"  this  electrode  arrangement  was  adopted  as  an 
approach  to  the  problem  of  a  single  sphere  in  an  infinite 
medium.  The  assumptions  made  have  been  discussed  in  the 
"Introduction" . 

A.  Electric  Field 

1.2  General  Solution 


The  electric  field  satisfies  Laplace's  Equation,  as 
shown  in  Appendix  A  , 

V2U  =  0  1.2.1 


subject  to  the  boundary  conditions 


and 


U  (a)  =  UQ 

U (b)  =  0 


1.2.2 

1.2.3 


The  Laplace  operator  in  spherical  co-ordinates  is 

9  1  a  9  a  la  a2 

v  =  — 7c  — —  (r  — —)  +  -iy  —(sin  0  — ) 

r  9r  ar  r  sin  636,  96 
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“2 — - 

r  sin^  6  a  4> 


, 

airf3  ‘  aolcfrowbpi-  n  I  9  . 


is  io  maldon 


xion^q-qA  nl  woria 
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16  ? 
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where  r  is  the  radius,  <j>  the  longitude  and  0  the  colatitude. 

From  symmetry  considerations  there  is  no  6  or  <f>  de¬ 
pendence  in  this  problem.  Therefore,  equation  1.2.1  becomes 

13  7  3U 

—  (r^  — )  =  0 


3  r 


3  r 


Multiplying  by  r  and  integrating  yields 


3  U 


3  r 


=  C 


1 


Dividing  by  r  and  integrating  again  gives 

C. 


U  =  - 


'1 


+  C. 


1.2.4 


This  equation  is  the  general  solution  for  U(r). 


1 . 3  Application  of  Boundary  Conditions 


The  constants  and  are  determined  by  applying 
the  boundary  conditions.  From  equations  1.2.3  and  1.2.4  it 
follows  that 

C1 

0  = - -  +  c9 

b 


or 


From  equations 


1.2.2  and  1.2.4 


1.3.1 
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C1  = 

u 

o 

1  _  1 
b  a 

II 

U  ab 
o 

a-b 

Inserting  these  values  in  equation  1.2.4  yields 


U  = 

U  ab  1  U  a 

(  °  )  °  1.3.2 

b-a  r  b-a 

Equation  1.3.2  describes 

the  potential  distribution  for  the 

problem. 

1.4  Evaluation  of  the  Heat  Generation  Term  in  the 
Diffusion  Equation 

As  shown  in  Appendix  A,  the  heat  generation  term  in 


the  diffusion  equation  is 

given  by  the  expression 

E  • 

l 

1.4.1 

The  electric  field  vector  may  be  expressed  as 

-> 


E  = 

=  -vu 

Considering  symmetry , this 

yields 

E 

r 

3  U 

3r 

Applying  this  to  equation  1.3.2  yields 
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Ohm's  Law  in  differential  form  is 


1  =  a  E 

e 


The  heat  generation  term  (1.4.1)  can  therefore  be  stated  as 


E  •  l 


CT 


-S.  1  E  I  2 


O 


Substituting  for  |e|  from  1.4.2  gives 


->  •> 

E  •  i 


a  U  ab  0  1 

_  =  _£  (_° _ )2  _ 

,  1  4 

b-a  r 


For  convenience  define  a  parameter  G  by  the  relation 


a  U  ab  „ 
G  =  —  (^)2 

°t  b_a 


With  this  definition  the  generation  term  in  the  diffusion 
equation  is  given  by 


1 


1.4.3 


B.  Thermal  Field 
1 . 5  Separation  of  Variables 

The  thermal  field  satisfies  the  diffusion  equation 

with  a  heat  generation  term  added,  as  shown  in  Appendix  A. 

In  mathematical  form 

~  E  •  i  1  3T 
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Using  1.4.3  this  becomes 

9  G  1  3T 

V  T  +  -j.  = - 

r  k  3 1 


Once  again  symmetry  rules  out  the  6  and  cf>  dependence.  The 
basic  partial  differential  equation  for  T(r,t)  is  therefore 


8  T 


1  3  T 


r2  3r  3 r  r4  k  3t 


1.5.1 


For  the  purpose  of  solving  this  equation  the  method 
of  separation  of  variables  is  employed.  Assume 


T(r,t)  =  R (r)  S(t)  +  f(r)  1.5.2 

where  R(r)  and  f(r)  are  dependent  only  on  radius  and  S(t)  is 
a  function  of  time  alone.  Equation  1.5.1  then  becomes 
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3  r 

(r  — ) 

3  r 
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3r  r 
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The  function  f  may  now  be  chosen  to  satisfy 


1  3 


3  f 
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r  3  r  3  r 
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and  we  are  left  with 


S  3 

—  —  (r 
r  3  r 


3  R 


R  3  S 


3  r  k  3  t 


which,  when  divided  by  RS ,  becomes 
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Since  the  left  hand  side  of  this  equation  is  a  function  of  r 
alone  and  the  right  hand  side  a  function  of  t  alone,  the  most 
general  condition  possible  is 


1 

9 

0  9  R 

1 

9  S  n 

o 

— 

(r2  — ) 

=  — 

—  =  ±A2 

1.5.5 

Rr 

9  r 

9  r 

kS 

9  t 

where  X  is  a  real  constant. 


1 . 6  Steady  State  Solution 


The  function  f  describes  the  final  or  steady  state 
temperature  distribution  in  the  medium  between  the  spheres. 
This  follows  from  1.5.2  since  S(t)  must  vanish  as  t  approaches 
infinity . 

Consider  equation  1.5.3 


1  9  9  9  f 

-o  —  (r  — )  + 
r  9  r 


=  0 


(1.5.3) 


9r  r 

2 

Multiplying  by  r  and  integrating  yields 

-  9f  G 

r  —  =  -  +  C. 

9r  r 

2 

Dividing  by  r  and  integrating  once  more,  the  general  solu¬ 
tion  for  f  is  obtained. 
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1 . 7  Time  Dependent  Solution 


Consider  the  time  dependent  equation  contained  in 
1.5.5,  viz . 

1  3S  2 


=  ±A' 


kS  3 1 


Rearranging  the  terms  yields 

3  S  ~ 

-  =  ±A  k3t 

S 

which,  when  integrated,  results  in  the  solution  for  S 


In  S  =  ±X  kt  +  In  C, 


The  nature  of  the  problem  requires  that  the  solution  remain 
bounded  as  time  tends  to  infinity,  forcing  selection  of  the 
negative  sign.  This  gives  as  the  solution  for  S(t) 


S(t)  =  c5  e 


-  A  2kt 


1.7.1 


1 . 8  Radius  Dependent  Solution 


Selecting  the  radius  dependent  equation  from  1.5.5 

13  ~  3  R  ~ 

- 2  —  (r  — )  = 

Rr  3  r  3  r 

On  multiplying  by  R  and  collecting  all  terms  on  the  left  this 
becomes 

3  2R  2  3R  „ 

- rr  +  —  —  +  A  R  =  0  1.8.1 

3  r  r  3r 

which  may  be  transformed  into  Bessel's  Equation  by  the  follow- 
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ing  procedure 
Define 


Q  =  r  p  R 


1.8.2 


and 


y  =  cr 


1.8.3 


The  "inverse"  relations  are  then  given  by 


R  =  r^  Q 


1.8.4 


and 


,V/z 

r  =  (-)  ' 

c 


1.8.5 


Differentiating  1.8.4  with  respect  to  r  yields 


9  R  ,  9Q 

—  =  prP'h  +  rP  — 


9r 
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and 
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3  r  9r 


Substituting  in  1.8.1  for  R  and  its  derivatives  gives 
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By  the  chain  rule  of  differentiation,  using  1.8.3 
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Similarly , 
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3Q 
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Again  using  1.8.3 
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Inserting  these  relations  in 

2 

1.8 

. 6  results  in 
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2~) 
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Multiplying  by  r 

3  2Q 
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2  2  2z  7 

c  z  r  — j  +  (2p+2)czr  — 

9  y  9y 


the  above  equation  becomes 
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p(p-l)  +  2p  +  X^r^ 
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0 


Substitute  in  this  relation  the  expression  for  r  given  by 
lo8.5.  This  yields 
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2  2  y  ^  Q  y  9  q 

c  z  —  2  +  (2p+2)  cz  —  —  + 

c  dy  c  3y 


p(p-l)  +  2p  +  X' 
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Q  =  0 
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Q  =  0 

1.8.7 


The  desired  form  of  equation  1.8.7  is 

BV 
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9  a2v  .v  9  9 

— +  x  —  +  (x  -  n  )v 
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ax' 


ax 


This  form  is  obtained  by  choosing  z  =  1,  p  =  ,  and  c  =  X. 

The  relation  1.8.7  thus  becomes 
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2  /  X  2 

y  -  (-) 

2  - 


Q  =  0 


3y  3y 

which  is  Bessel's  Equation  of  order  h.  The  general  solution 
of  this  equation  is 


Q  =  C6J,  (y)  +  C7J  ,  (y) 


1.8.8 


where  0/  (y)  and  J  ±  (y)  are  the  Bessel  functions  of  order  h. 
^2  - $ 

Combining  1.8.2,  1.8.3  and  1.8.8  gives 
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Applied  to  equation  1.8.9  these  result  in  the  following  as 
general  solution  for  R(r): 


R  (r) 


=  /2_ 

/  TT  X 


sin  Xr 


cos  Xr 


+  C. 


1.8.10 


Having  obtained  the  solutions  for  each  of  the  components 
of  equation  1.5.2  it  is  now  possible  to  write  the  general  so¬ 
lution  for  T(r,t).  The  original  equation  (1.5.1)  has  an  in¬ 
finite  number  of  solutions  corresponding  to  various  values  of 
the  separation  constant  X.  The  general  solution  thus  involves 
a  Fourier  series,  viz.  that  portion  of  the  following  expres¬ 
sion  contained  within  the  summation.  Inserting  the  solutions 
1.8.10,  1.7.1,  and  1.6.1  into  1.5.2  gives,  as  the  general 


solution  for  T 

T  (r  ,  t)  = 


kt 
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Cc  sin  (X  r) 
6n  v  n  ' 


+  C_  cos  (X  r) 
7n  n 


TT - +  C4 

2r  r 


The  arbitrary  constants  may  be  incorporated  into 

the  constants  Cr  and  C_  with  no  loss  of  generality.  This 
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It  now  remains  to  choose  the  unknown  constants  A  , 

n 

B  ,  and  A  ,  in  such  a  fashion  as  to  satisfy  the  initial  and 
n  n 

boundary  conditions. 


1 . 9  Application  of  the  Boundary  Conditions  on  T 


The  following  requirements  are  to  be  imposed  on  the 
solution  1.8.11. 


1) 
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2)  T  (b ,  t )  =  T 


amb 


Applying  condition  1  to  1.8.11 
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The  restriction  "for  all  times"  implies  that  the  time  in¬ 
dependent  and  time  dependent  portions  must  each  be  zero. 
That  is 


3  +  2  0 

a  a 


or 
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and 
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n=0 


2 


TT  A 
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-A 
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kt 


A  sin  A  a  B  cos  A  a 
n  n  n  n 


+ 


n 


(A  cos  A  a  -  B  sin  A  a) 
n  n  n  n 


=  0 


The  summation  will  certainly  vanish  if  each  of  its  terms 
vanishes.  Therefore,  it  is  required  that 


1 


A 


n 


-  — (A  sm  A  a  +  B  cos  A  a)  +  —  (A  cos  A  a 
2  n  n  n  n  n  n 

a  a 


-  B  sin  A  a)  =  0 

n  n 


for  all  n. 


Rearranging  and  multiplying  by  a  yields 


A  (A  a  cos  A  a  -  sin  A  a)  B  (A  a  sin  A  a 

n  n  n  n  n  n  n 


+  cos  A  a)  =  0 
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Applying  condition  2  to  1.8.11  gives 
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(A  sin  A  b  +  B  cos  A  b) 
n  n  n  n 


- 7  +  —  +  Cd 

2b  ab 


The  transient  portion  must  be  zero  because  the  left  hand  side 

has  no  time  dependence.  It  therefore  follows  that 
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and 


-  Xnkt 


^  ^  n\i  b 


n=0 


(A  sin  X  b  +  B  cos  X  b)  =  0 

n  n  n  n 


Again  each  term  of  the  summation  must  be  zero.  That  is 


A  sin  X  b  +  B  cos  X  b  =  0  for  all  n  1.9.2 

n  n  n  n 


We  now  have  applied  two  conditions  and  therefore  ex¬ 
pect  to  have  lost  two  degrees  of  freedom.  The  initial  con¬ 
dition  yet  to  be  applied  requires  that  one  degree  of  freedom 
be  retained  in  the  amplitudes  A  and  B  .  This  can  be  accom- 
plished  only  if  the  two  equations,  1.9.1  and  1.9.2,  are 
linearly  dependent.  We  therefore  specify  that  the  determin¬ 
ant  of  the  coefficients  of  A  and  B  be.  zero.  In  matrix  form 

n  n 

these  equations  are  written 


— * 

—  — 

— 

X  a  cos  X  a  -  sin  X  a 

-X  a  sin  X  a  -  cos  X  a 

A 

0 

n  n  n 

n  n  n 

n 

sin  X  b 

cos  X  b 

B 

0 

n 

n 

n 

Setting  the  determinant  to  zero  gives  the  following  equation. 

X  a(cos  X  b  cos  X  a  +  sin  X  b  sin  X  a) 
n  n  n  n  n 

t(sinXbcosXa-cosXbsinXa)  =  0  1.9.3 

n  n  n  n 

Using  the  trigonometric  identities 

cos  (x-y)  =  cos  x  cos  y  +  sin  x  sin  y 
sin(x-y)  =  sin  x  cos  y  - 


cos  x  sm  y 
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equation  1.9.3  may  be  written 


X  a  cos  X  (b-a)  +  sin  X  (b-a)  =  0  1.9.4 

n  n  n 

The  eigenvalues  are  defined  by  equation  1.9.4.  For 
the  purpose  of  numerical  evaluation  it  is  convenient  to  in¬ 
troduce  the  variables  a  ,  defined  by  the  relation 

n '  2 


That  is 


a 


n 


X  (b-a) 
n 


b-a 


The  variables  may  be  regarded  as  the  eigenvalues  in  dimen¬ 
sionless  form. 

Substituting  for  Xn  in  1.9.4  and  dividing  by  cos 


yields 


tan  a 

n 


aa 

n 

b-a 


1.9.5 


This  equation  is  an  alternate  form  of  1.9.4,  the  independent 
variable  in  this  case  being  dimensionless.  Figure  1.2  il¬ 
lustrates  the  positions  of  the  roots  of  1.9.5. 

It  should  be  noted  that,  since  X  =  0 ,  we  must  set 

'  o 

A  =  B  =  0  to  avoid  indeterminacy  in  T. 
o  o 

The  first  degree  of  freedom  has  been  lost  in  defin¬ 
ing  the  eigenvalues.  The  second  has  been  lost  in  defining 
the  ratio  B^/A^ ,  leaving  the  third  in  the  amplitudes  as 
required . 

For  the  purpose  of  applying  the  initial  conditions 
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Positions  of  the  Roots  of  the  Eigenvalue  Equation , 
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Argument  a  in  Radians 


1.17 


it  is  necessary  to  find  a  function,  say  ,  which  is  ortho¬ 
gonal  over  the  interval  (a,b)  and  satisfies  both  boundary 
conditions.  The  orthogonality  requirement  is  imposed  be¬ 
cause  we  wish  to  expand  an  arbitrary  function,  i.e.  the  ini¬ 
tial  condition,  which  must  also  satisfy  the  boundary  condi¬ 
tions,  in  a  Fourier  series.  Orthogonality  facilitates  deter¬ 
mination  of  the  coefficients  of  the  series. 

Consider  the  following  definition  of  W^: 

C  (sin  X  b  cos  X  r  -  cos  X  b  sin  X  r) 

TT  n  n  n  n  n 


where  Cn  is  an  arbitrary  constant.  Using  the  trigonometric 
identities  stated  above  this  equation  may  be  put  into  the 
more  compact  form 

C 

W  =  —  sin  X  (b-r) 

n  n 

r 

Thus  defined,  satisfies  the  boundary  conditions;  that  is 
for  r  =  b, 

C 

W(b)  =  —  sin  0 
b 

=  0 

and,  for  r  =  a, 


3W 

n 

C 

n 

aX  cos  X  (b-a)  +  sin  X  (b-a) 

3r 

a^ 

r=a 

n  n  n 

=  0 

according  to  the  eigenvalue  equation,  1.9.4. 
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This  function  is  also  orthogonal  over  the  interval 
(a,b),as  shown  in  Appendix  B. 

The  equivalence  of  and  the  radius  dependent  por¬ 
tion  of  the  n^*1  term  of  the  general  solution  for  T  may  be 
verified  by  choosing 


n 


■C  cos  X  b 
n  n 


X  TT 
n 


and 


B 


n 


C  sin  X 
n 


b  / 
n  / 


/x  rT 

n 


The  following  solution  for  T(r,t)  has  now  been  es¬ 
tablished  . 


T  ( r ,  t)  = 


~  -XR  kt  C 

}  e  —  sin  X  (b-r) 

n=l  r 


- +  -  +  T  ,  +  - 7T  —  - 

„  2  amb  2  , 

2r  ar  2b  ab 


1.9.6 


1.10  Application  of  the  Initial  Condition  on  T 


The  constants  Cn  are  still  unknown  and  may  be  deter¬ 
mined  by  applying  the  initial  condition.*  We  require 


T (r , 0 )  =  T 


amb 


From  1.9.6 


it  is  at  this  point  that  other  initial  conditions  may 
be  applied.  See  section  5.2. 
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sin  A  (b-r) 
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2r  ‘ 
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or 
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sin  A  (b-r)  =  — 
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n=l  r  2r 
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Multiply  both  sides  by  r  sin  A  (b-r)  and  integrate  from  a 
to  b .  By  the  orthogonality  property  (Appendix  B)  this  re¬ 
sults  in 


C  7  b  G 

—  (b  -  a  sin  A  (b-a) )  =  /  r  sin  A  (b-r)  ( - ~ 

n  n  J  n  '2 

2  a  2r 
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Solving  for  C 
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b-a  sin  A  (b-a)  a 
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,b  G 

/  r  sin  An  (b-r)  (■ 


2r‘ 
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1.10.1 


Therefore  the  general  solution  is  given  as 


T(r,t)  = 


-A  kt 
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y 

n=l 


C  sin  A  (b-r) 
n  n 
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1.11  Conclusion 

One  of  the  integrals  involved  in  determining  the  co- 
b  sin  X  (b-r) 

efficients,  viz.  /  - - -  dr,  is  transcendental.  For 

a  2r 

this  reason,  and  because  the  eigenvalue  equation  requires 
iterative  techniques  for  its  solution,  the  use  of  a  digital 
computer  becomes  virtually  essential.  Chapter  4  is  concerned 
with  the  computational  problems  involved  in  such  a  solution 
and  presents  results  which  were  obtained  using  a  program 
given  in  Appendix  D. 

The  procedure  beginning  with  equation  1.8.2  of  trans¬ 
forming  the  radius  dependent  partial  differential  equation 
into  the  standard  form  of  Bessel's  equation  can  be  avoided 

g 

by  a  technique  suggested  in  Carslaw  and  Jaeger  .  This  al¬ 
ternate  approach  is  used  in  Chapter  2  and  ultimately  results 
in  the  same  form  of  solution  for  the  bounded  problem.  Both 
methods  have  been  employed  primarily  to  illustrate  the  al¬ 
ternatives  available. 
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CHAPTER  2:  SINGLE  SPHERE  IN  AN  INFINITE  MEDIUM 


2.1  Introduction 


In  Chapter  1  the  solution  for  the  electric  and  thermal 
fields  between  two  concentric  spherical  shells  was  developed. 
If  the  radius  of  the  outer  shell  is  allowed  to  increase  with¬ 
out  limit  it  becomes  necessary  to  adopt  a  somewhat  different 
technique.  Chapter  2  deals  with  this  case;  the  analytic  de¬ 
scription  of  the  electric  and  thermal  fields  surrounding  a 
single  spherical  shell  embedded  in  an  infinite  medium.  Ini¬ 
tial  and  boundary  conditions  and  the  assumptions  made  in  the 
solution  are  discussed  in  the  "Introduction"  (page  ii) . 

A.  Electric  Field 

2 . 2  Application  of  Boundary  Conditions 


The  general  solution  for  the  steady  state  electrical 
potential  for  a  spherically  symmetric  electrode  arrangement 
was  developed  in  Chapter  1  and  is  given  by  equation  1.2.4 

^1 

U  - - -  +  c2  (1.2.4) 

The  boundary  conditions  for  the  present  problem  are 


U (a)  =  Uq  2.2.1 
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Applying  2.2.2  to  1.2.4  results  in 


lim  U 

2f->-oo 


lim  (-  —  +  C2) 

r->°°  r 


-  C2 
=  0 

That  is  C2  =  0 

Applying  2.2.1  yields 

ci 

u  - - - 


or 


C1  =  -aUo 

The  potential  distribution  is  therefore  described  by 

U  a 

U  (r)  =  —  2.2.3 

r 


2.3  Evaluation  of  the  Heat  Generation  Term  in  the 
Diffusion  Equation 


It  was  shown  in  Section  1.4  that  the  required  term  may 
be  written 
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and  that 
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Applying  this  to  equation  2.2.3  yields 
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which  gives  as  the  heat  generation  term 

(E • 1 )  a  U *  2a2 

e  o 

4 

at  at  r 

For  convenience  a  parameter  G,  defined  by  the  following  re¬ 
lation,  is  introduced: 


Note  that  this  definition  of  G  is  consistent  with  that  in  sec¬ 
tion  1.4  when  the  outer  radius  b  is  allowed  to  increase  with¬ 
out  limit. 

The  required  heat  generation  term  in  the  diffusion 
equation  thus  becomes 


(E-i)  G 


2.3.1 


B.  Thermal  Field 
2 . 4  Separation  of  Variables 

The  following  equation  describes  the  thermal  field 

~  E • 1  1  9T 

V  T  +  -  =  -  - 

o  ^  k  9 1 

Using  2.3.1  the  partial  differential  equation  governing  T 
becomes 
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Assuming  radial  symmetry  this  relation  reduces  to 
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d2T  2  3T  G  1  3  T 

- 2  +  —  —  +  — -r  —  —  —  2.4.1 
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Following  the  suggestion  in  Carslaw  and  Jaeger  a  function 
W(r,t)  is  defined  by  the  following  relation 

W(r,t) 

T  =  -  2.4.2 
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Then 
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Substituting 
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relations 


2  3W  1  32W 

~~2  +  —  o’ 
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in  2.4.1  and  multiplying  r  yields 
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2.4.3 


Employing  the  method  of  separation  of  variables,  the  functions 
f(r)  and  V(r,t)  are  defined  such  that 


W(r,t)  =  V(r,t)  +  f(r)  2.4.4 

^In  terms  of  V  and  f  equation  2.4.3  is 
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Since  V(r,t)  is  time  dependent  and  the  general  solution 
must  remain  bounded  as  time  tends  to  infinity,  the  function  V 
must  ultimately  vanish.  This  implies,  considering  the  above 
equation,  that 
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and  that 
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2.4.6 


The  solution  of  these  two  equations,  subject  to  the  boundary 
conditions,  will  now  be  effected. 

\ 

2 . 5  Steady  State  Solution 

The  steady  state  distribution  VJ(r,x)  is  described 
by  the  function  f,  which  satisfies  2.4.5 
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Integrating  twice  gives 
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This  equation  is  the  general  solution  of  the  steady  state  dis¬ 
tribution  , 

f  =  W  ( r ,  °° ) 
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2.6  Derivation  of  the  Modified  Boundary  Conditions, 
Application  to  the  Steady  State  Solution 


There  are  two  boundary  conditions  on  T 
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2)  lim  T(r,t) 
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3T 

3 r  r“  r  3r 

or,  in  terms  of  V  and  f,  using  2.4.4 
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2.6.2 


The  function  V  is  time  dependent  whereas  f  is  not. 

It  therefore  follows  that  each  function  must  satisfy  the 
boundary  conditions  independently.  Caution  must  be  exercised 
however,  to  avoid  changing  the  sign  of  one  part  without  chang 
ing  that  of  the  other.  Mathematically  this  is  expressed  as 


V  13V 
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To  determine  the  constant  C ^  in  the  solution  for  f 
the  condition  2.6.4  is  applied.  Recalling  the  general  solu¬ 
tion  for  f,  2.5.1 

G 

f  =  +  c]_r  +  c2  (2.5.1) 

2r 


then 

3  f 
3r 

Substituting  in  2.6.4  gives 
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a 

Boundary  condition  2  may  be  used  to  determine  C-^. 
From  2.4.2  and  2.4.4  it  follows  that 

V  f 


which,  by  the 
two  equations 
is 
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time  dependence  argument,  can  be  separated  into 
for  application  of  boundary  condition  2.  That 
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Inserting  the  function  f  we  obtain 


lim  (-  — 
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That  is 
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The  complete  solution  for  f  is  therefore 
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Consider  the  transient  portion  of  the  function  W. 

For  the  purpose  of  determining  the  function  V(r,t) 
it  is  convenient  to  introduce  a  function  S(r,t)  such  that  the 
inner  boundary  condition  on  S  is  homogeneous.  Equation  2.6.3 
states  that 
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If  S  is  defined  by  the  relation 
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the  inner  boundary  condition  on  S  is 

S (a , t)  =  0 

which  is  homogeneous. 

The  differential  equation  governing  S(r,t)  may  be 


found  as  follows. 
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Differentiating  equation  2.6.6  yields 
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Differentiating  once  more 
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Comparing  this  with  2.6.7  shows  that  S  satisfies  the  same 
differential  equation  as  V,  namely 
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2 . 7  Application  of  the  Initial  Condition 
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or,  in  terms  of  V  and  f. 
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Multiplying  by  r  and  rearranging, 
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From  2.7.1  it  follows  that 
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Recall  the  definition  of  S 
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Therefore,  using  2.7.1  and  2.7.2 
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For  convenience  we  define 


x  =  r  -  a 
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or 


r  =  x  +  a 


r  >  a 


Let 


Q (x,t) 


S  (x  +  a  ,  t) 


V 


V  ip  amxPd  ni 

' 


dims 


T* 


^pnipnexTSSi  fcftfi  d  i  v  i- •  T ^  X  c/M .. 


-  J 


Tx  *'  - 


dm* 


xs; 


i  / 


dBiid  ewolloi  di  I.  fno'1^ 


—  rt 


0»d 


* 


/•* 

’ 


2  "io  xioxixnxiab  »rid  IIsppH 


X  s 


V 


S.r/S.Ms  X.t'/S  pniao  .sxoisxPfiT 


"x£ 


a  -  x 


l 


i 

6  JS  1 


a  +  x  *  x 


deJ 


■ 


2.11 


Then,  at  time  t  =  0 

-G  G  aG 

Q.  ( x  ,  0 )  =  -  +  —  -  - j  x  0  2.1.  A 

2  (x+a)  a  2  (x+a) 

The  problem  now  is  to  find  Q(x,t)  which  must  satisfy 
the  diffusion  equation  and  the  initial  condition  2.7.4. 

As  shown  in  Appendix  A  the  function 


A 


-(x-x') 

4kt 


2  /ktir 

is  a  particular  solution  of  the  diffusion  equation  in  x'  and 
t.  "A"  is  an  arbitrary  constant.  Since  the  equation  is 
linear,  any  sum  of  particular  solutions  is  also  a  solution. 
Hence,  if  we  regard  A  as  a  function  of  the  dummy  variable  x' 
and  define  A(x')  such  that  it  is  identical  to  Q(x',0),  the 
solution  for  Q(x,t)  may  be  expressed  as  follows: 
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It  can  be  shown  that  the  limit  of  the  right  hand  side 
of  equation  2.7.5  as  t  approaches  zero  is  Q(x,0). 

In  order  to  assure  convergence  of  the  integral  in 
2.7.5,  we  define  Q  in  the  interval  x^  0  such  that 


Q  (x  <  0 , 0)  =  -Q  (x  >  0 , 0) 


Explicitly,  Q  is  then  given  by 


II  .  s 


Dis 


3- 


.  K. .  S 


e  \{&+x)  £ 


v 


noiddmjl- arid  A  xibneqq A  ns  nwoi'a 


A { ’ x-  x)  - 


A 

0  - 


wJSNs 


noiifix^  art;  ©oniB  .^ns^anoo  y*®'*^*6  ne  ei  A  #i 


.aoiJWos  2  osIb  ai  enoi^fcite  MXuoiiieq  5o,mo8  <f<is  tM»£U 

B 

‘  A 

jaw9riol  M.biaBSiqpt*  ad  Yem  noxds/Ioit 


«>  + 


(d,x)0 


' 


x  L&ivosttX  sfli  Xo  eonogievnoo  smsifi  02  n.i 


■  C  ,0  -  :  , 0-  -  '3 * ;  -  /0 


Y<J  n-evip  xts/i^  ai  2  .  qx3 


2.12 
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Q(x,0) 
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Using  these  equations  in  2.7.5  yields 
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Consider  the  first  of  these  integrals. 

Let  y  =  -x ' 

then  dy  =  -dx ' 

and  the  integral  becomes 
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Here  y  is  a  "dummy"  variable.  A  comparison  of  this  with  the 
second  integral  reveals  that 
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Q(x,t)  = 
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Evaluation  of  this  integral  must  be  done  numerically. 

2 . 8  Inversion  of  Transformations 

Having  obtained  the  function  Q  we  now  work  back 
through  the  various  transformations  to  find  T(r,t). 

Q(x,t)  is  related  to  S(r,t)  by  the  following  change 
of  variable: 


That  is 


x  =  r  -  a 


S  (r,  t)  =  Q  (r  -  a,  t) 


The  function  S  may  therefore  be  expressed  as 
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S  is  related  to  V  by  the  equation 


S  =  -V  +  a 
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or,  dividing  by  +a, 


V  3V 
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The  homogeneous  solution  of  this  equation  is 


V  =  C  e 


r/a 


Since  V  must  remain  bounded  as  r  tends  to  infinity  (boundary 
condition  2) ,  it  follows  that 

C  =  0 

The  particular  solution  is  obtained  by  the  following 

means : 
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The  right  hand  side  is  the  derivative  of 
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Therefore,  taking  the  definite  integral  of  both  sides  between 
the  yet  undefined  limits  g  and  h 


_r/a  T  T 

e  '  V 


=  f 

J 

g  g 


h  e-r/a  S 


dr 


2.8.3 


.«+  vd  pnjLbivi:& 


Vf  V 


4 


'  J. 


■£'  D  a  V 


■ 

dsdd  yjwolloS  dX  -■»  (&  .  next ~  tbflO'-' 

D-:  : 


gniwoliol  ai-d  .;3  banifiddb  li  t  roiit  loa  is  I  • 


:  ensainu 


S ,  B .  £  noXdsnpa  y  IqidluM 

>*.  i‘ 


,v..  i  -  l 

*VfS  .  5  /  9 


•s 


s\l~ 


©8 


45 


it  is  '  t  ■  *.  '  ^fcri  ,:s  : 


lS>c>wdad  gabis  ridod  'io  ®iXaX$ab  arid  Pdi>6"  ,®*ie-t9*ionT 


, 


rf  bns  9  adinul  banidabnu  day  add 


2.15 


To  determine  the  limits  g  and  h  we  use  the  remaining 
boundary  condition,  viz. 

lim  T  =  T 


2T->°o 


amb 


2.8.4 

(2.6.2) 


and  the  fact  that  r  should  be  retained  as  the  independent  space 
variable.  The  former  restriction,  however,  is  imposed  on  the 
function  T  whereas  the  limits  in  2.8.3  are  involved  in  evalu¬ 
ating  V(r,t) .  To  determine  g  and  h  it  is  therefore  necessary 
to  find  the  condition  on  V(r,t)  which  is  equivalent  to  2.8.4. 
This  is  accomplished  as  follows: 

From  2.4.2  and  2.4.4 
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or,  rearranging  and  multiplying  by  r 
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equation  2.8.3  will  yield 
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The  function  W  was  related  to  V  by 
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The  original  transformation  was 
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The  complete  solution  for  T  is  thus 
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where  S  is  given  by  the  equation  2.8.1. 


2 . 9  Conclusion 

In  principle  the  analytic  description  of  the  thermal 
field  is  given  by  2.8.6.  This  equation  involves,  however,  a 
double  integration  which  must  be  done  numerically,  by  means 
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of  a  digital  computer.  The  difficulties  involved  in  such  a 
solution  are  discussed  and  results  based  on  a  program  given 
in  Appendix  E  are  presented  in  Chapter  4 . 


3.1 


CHAPTER  3 .  CONCENTRIC  CYLINDERS 

3 . 1  Introduction 

This  chapter  deals  with  the  analytic  description  of 
the  electric  and  thermal  fields  between  two  infinitely  long 
co-axial  cylinders.  The  development  parallels  that  of  Chap¬ 
ter  1  and  incorporates  essentially  the  same  assumptions  and 
conditions . 


A.  Electric  Field 


3 . 2  General  Solution 


Laplace's  equation  in  cylindrical  coordinates  governs 
the  electric  potential: 
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Symmetry  dictates  no  0  or  z  dependence,  hence  the  above  equa¬ 
tion  reduces  to 
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Rearranging  and  integrating  gives 
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s  equivalent  to 


U  =  C1  In  r  +  C2  3.2.2 

Equation  3.2.2  is  the  general  solution  for  the  electric  po¬ 
tential  distribution. 


3 . 3  Application  of  the  Boundary  Conditions 
The  conditions  to  be  applied  are 
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Inserting  3.3.1  into  3.2.2  gives 
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In  a  -  In  b 


or,  explicitly 
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In  a  -  In  b 


The  particular  solution  for  U(r)  is  therefore 


(In  r  -  In  b) 
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3.3.3 


°  (In  a  -  In  b) 


3.4  Evalution  of  the  Heat  Generation  Term  in  the 
Diffusion  Equation 

The  heat  generation  term  in  the  diffusion  equation, 
as  shown  in  Appendix  A,  is  given  by 


Using  Ohm's  Law  and  taking  account  of  symmetry  this  becomes 
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The  electric  field  intensity  is  given  by  the  expression 

9  U 


Using  3.3.3  the  generation  term  can  therefore  be  written 
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3.4 


For  convenience  a  parameter  G  is  defined  as  follows 


G 


a  U  ‘ 
e  o 


at  In' 


a 
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3.4.1 


With  this  definition  the  generation  term  in  the  diffusion  equa¬ 
tion  is  given  by 

G 

~2 
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B.  Thermal  Field 
3 . 5  Separation  of  Variables 


Under  the  present  conditions  of  cylindrical  symmetry 

the  diffusion  equation  describing  the  thermal  field  is 

3* 2 3 * ST  1  3T  G  1  3T 

- +  —  — —  +  — j  =  —  —  3.5.1 

3r  r  3r  r  k  3t 

The  solution  of  this  equation  may  be  effected  by  the 
method  of  separation  of  variables.  Assume 


T  ( r ,  t )  =  R(r)  S(t)  +f(r)  3.5.2 

where  R(r)  and  f(r)  are  functions  of  radius  only  and  S(t)  is 
solely  a  function  of  time.  Equation  3.5.1  then  becomes 

3  2R  S  3  R  3 2  f  1  3  f  G  R  3  S 

S  - ^  t  —  -  - v  —  —  t  — 2  =  —  -  3.5.3 

3r  r  3 r  3r  r  3r  r  k  3t 


The  solution  for  T(r,t)  must  remain  bounded  as  time 
increases  without  limit.  For  this  reason  the  function  S(t) 
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3.5 


must  tend  to  zero.  Examination  of  the  above  equation,  3.5.3 
therefore  reveals  that  f,  the  steady  state  temperature  dis¬ 
tribution,  must  satisfy 

3 2  f  1  3f  G 

- 2  +  —  —  +  — 2  ~  0  3.5.4 

3r  r  3r  r 

Equation  3.5.3  then  becomes 

32R  S  3R  R  3S 

3r2  r  3r  k  3t 

Dividing  by  RS  and  rearranging  puts  this  into  the  form 

1  3  2R  1  3R  1  as 

R  3 r 2  Rr  3r  Sk  at 

Since  the  two  sides  of  this  equation  are  functions  of  two  in 
dependent  variables  the  most  general  case  possible  is  that 
each  side  is  equal  to  a  constant.  That  is 

1  3  2R  1  3R  1  3S  9 

- 9  + - = - =  ±  A  3.5.5 

R  3r  Rr  3r  Sk  3t 

The  problem  has  been  reduced  to  solving  equations 
3.5.4  and  3.5.5.  The  details  of  the  solutions  are  presented 
in  the  following  sections. 

3 . 6  Steady  State  Solution 

Consider  first  equation  3.5.4 
32f  1  3f  G 

— -rr  +  —  —  +  — —  0  (3.5.4) 
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Multiplying  by  r 
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Let 

or 

Then 
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Substituting  in 

That  is 


yields 
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3.6.1  gives 
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Integrating  twice  yields 

9  f 

=  -  Gy  +  C 

9y 


f  =  -  +  C,y  +  C? 

2  1  Z 

Substituting  for  y  from  3.6.2  results  in  the  following  general 
steady  state  solution. 

G  In2 * * S  r 

f  =  -  +  C1  In  r  +  C9  3.6.3 

2 *  1  ^ 


The  constants  and  may  be  determined  by  applying  the 
boundary  conditions. 

3 . 7  Time  Dependent  Solution 


We  have,  from  equation  3.5.5 
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Rearranging  and  integrating 
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In  S  =  ± A 2  kt  +  In  C. 
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S  =  C3  e 


±X2kt 


The  S  function  must  be  bounded  as  time  tends  to  infinity, 
forcing  selection  of  the  negative  sign  in  the  exponent.  The 
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function  S  is  therefore  given  by 


S(t)  =  C3  e 


-X2kt 


3.7.1 


3 . 8  Radius  Dependent  Solution 


The  remaining  portion  of  3.5.5  is 
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1  d  R 
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R  3  r  ‘ 
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Rr  3r 


+  X  =  0 


3.8.1 


Multiplying  by  r  R  and  collecting  terms  results  in  the  equation 

3.8.2 
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If  we  now  let  Ar  =  y  and  hence 
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equation  3.8.2  becomes 

„  3  2R 
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3  R  „ 

y  —  +  y  R  =  0 

3  Y 


which  is  Bessel's  equation  of  order  zero.  The  general  solu¬ 
tion  of  this  equation  is 
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or,  in  terms  of  r 
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C„J  Ur)  +  CcY  (Xr) 
4  o  5  o 


3.8.3 


It  is  now  possible  to  write  the  general  solution  for 


T(r,t).  Since  equation  3.5.5  is  linear,  any  sum  of  solutions, 
corresponding  to  different  values  of  the  separation  constant  X , 
will  also  be  a  solution.  This  results  in  a  Fourier-Bessel 
series  involving  discrete  values  of  the  separation  constant. 
These  constants  or  eigenvalues  must  be  determined  by  applica¬ 
tion  of  the  boundary  conditions. 


The  amplitude  coefficients  and  C,_  may  be  com¬ 


bined  with  no  loss  of  generality.  There  will  be  an  infinite 
number  of  these  corresponding  to  the  various  eigenvalues. 

With  these  considerations  the  general  solution  for 
T(r,t)  becomes,  from  equations  3.5.2,  3.6.3,  3.7.1  and  3.8.3 
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kt 
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oo 
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T(r,t)  =  l  e 
n=0 
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In  r  +  C1  In  r  +  C2 
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Evaluation  of  the  eigenvalues  X  and  the  constants  A  and  B 

^  n  n  n 

is  based  on  the  initial  and  boundary  conditions. 

3 . 9  Application  of  the  Boundary  Conditions 


The  conditions  to  be  applied  are 
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and  2)  T(b,t)  =  T  , 

amb 

The  general  solution  for  T  consists  of  a  steady  state 
and  a  transient  portion.  The  latter,  being  time  dependent, 
must  eventually  decay  out  and,  in  general,  will  not  be  constant 
in  magnitude  at  a  given  radius.  In  order  to  satisfy  the  boun¬ 
dary  conditions  we  must  therefore  treat  the  two  portions  se¬ 
parately,  requiring 


8  (RS) 
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RS 


=  0 


3.9.1 
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r=b 


for  all  times 

0  3.9.2 


and  that 
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f (b)  =  T 
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3.9.3 


3.9.4 


The  total  expression,  T  =  RS  +  f ,  will  then  meet  all  the 
boundary  conditions  for  all  times.  Working  first  with  the 
steady  state,  from  3.9.3  and  3.6.3 
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Applying  3.-.  9.4  to  3.6.3  gives 


Gin  b 


+  G  In  a  In  b  +  C~  =  T  , 

2  amb 


That  is 


C  =  t 
^2  amb 


G  In  b 

+  -  -  G  In  a  In  b 


The  desired  solution  for  f  is  therefore 
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3.9.5 


Now  consider  the  transient  solution 
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Applying  3.9.1 
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non  non 


=  0 


This  will  be  satisfied  if  we  require 
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From  3.9.2 
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which  will  be  satisfied  if 


A  J  (X  b)  +  B  Y  (X  b)  =  0  3.9.7 

no  n  non 

The  situation  is  now  analogous  to  the  concentric  sphere  prob¬ 
lem  in  that  one  degree  of  freedom  must  be  left  in  the  ampli¬ 
tudes  A  and  B  in  order  to  satisfy  the  initial  condition. 
This  is  possible  if  3.9.6  and  3.9.7  are  linearly  dependent. 

Setting  the  determinant  of  the  coefficients  of  A  and  B  to 
r  n  n 

zero  gives  the  following  eigenvalue  equation. 


J ' ( X  a)  Y  ( X  b )  -  J  ( X  b)  Y ' ( X  a)  =  0  3.9.8 

on  on  on  on 


Define  a  function  W  (r)  as  follows 
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W  ( r )  =  C  J  '  ( X  a)  Y  ( X  r)  -  J  ( X  r)  Y'(X  a) 

n  ni^on  on  on  on_ 


3.9.9 


where  is  an  arbitrary  constant.  Thus  defined,  satis¬ 
fies  the  boundary  conditions  and  the  differential  equation. 
This  procedure  is  equivalent  to  choosing  the  constants  An  and 

B  as 
n 

A  =  -  C  Y' (X  a) 

n  non 


and 


B  =  C  J' (X  a) 
n  non 


The  solution  for  T  thus  becomes 


00  -X  ^kt 

T  (r ,  t)  =  l  e  n  W  (r)  +  f(r) 
n=0 


The  initial  condition  may  now  be  applied  to  define 
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3.10  Application  of  the  Initial  Condition 


The  initial  condition  specified  is 


T(r,0)  =  T 


amb 


Therefore,  employing  the  general  solution  for  T(r,t),  3.9.5 
and  3.9.9,  we  have 


1  ,  =  I  W  r  +  T 

amb  L n  n 
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or , rearranged 
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3.10.1 


The  orthogonal  properties  of  the  functions  are  employed 

to  evaluate  the  constants  C^.  As  shown  in  Appendix  C 
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Multiplying  both  sides  of  3.10.1  by  rW^  and  integrating  from 
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a  to  b  yields 


H  C  ‘ 
n  n 


That  is 
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=  /  rw 
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Gin  r 


-  G  In  a  In  r 


—In  b  +  G  In  a  In  b 
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1  b  rW 
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3.10.2 


Equation  3.10.2  determines  the  constants  C^.  The  complete 
solution  is  then 


°°  -X  kt  G  2 

T(r,t)  =  (  J  e  n  W_)  -  —  In  r  +  G  In  a  In  r 


n=0 


n 


+  —  In  b  -  G  In  a  In  b  +  T 


amb 


3.11  Conclusion 

The  solution  involves  transcendental  integrals  in 
the  equations  involved  in  defining  the  constants  C^.  The 
problem  was  programmed  but  time  limitations  did  not  permit 
inclusion  either  of  the  program  or  any  numerical  results. 

The  logarithmic  dependence  of  both  electric  potential 
and  temperature  rules  out  the  possibility  of  extending  the 
results  to  a  single  infinitely  long  cylinder  in  the  same  man¬ 
ner  as  the  concentric  spheres  problem  was  extended  to  that  of 
the  single  sphere.  Physically  either  of  the  cylinder  problems 
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is  impossible  since  both  require  an  infinite  quantity  of 
charge  to  realize.  The  results  are  of  interest,  however,  in 
that  the  fields  around  a  single  finite  cylinder  must  be  simi¬ 
lar  to  those  described  herein  under  certain  restrictions. 
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CHAPTER  4 .  NUMERICAL  TECHNIQUES 


4.1  Introduction 


This  chapter  is  concerned  with  the  numerical  tech¬ 
niques  and  problems  involved  in  programming  the  solutions 
given  in  Chapters  1  and  2  for  a  digital  computer.  The  numer¬ 
ical  analysis  methods  used  for  integrations  and  for  finding 
the  roots  of  the  eigenvalue  equations  are  discussed  along  with 
any  particular  difficulties  encountered  for  the  specific  prob¬ 
lems  in  question.  The  programs  themselves  will  be  found  in 
Appendix  D;  results  obtained  with  these  programs  are  tabulated 
in  Appendix  E  and  plotted  at  the  end  of  this  chapter. 

A.  Concentric  Spheres 

4.2  Simplification  of  Expression  Defining  Cn 


In  Chapter  1  the  theoretical  solution  for  the  concen¬ 
tric  sphere  problem  has  been  derived.  The  solution  found  re¬ 
quires  evaluation  of  the  constants  C^  which  are  defined  by 
the  following  equation. 
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2G 


b  -  a  sin  a 

n 


sin  A  (b-r) 
n 

2r 


dr 


b  sin  A  (b-r)  b  1 

-  /  - — -  dr  -  /  r  sin  A  (b-r)  dr  ( — ^ 

a  a  a  n  2b 

1 

ab 


(1.10.1) 
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Integrating  the  second  term  results  in 
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b  -  a  sin  a 
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b  sin  A  (b-r) 
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2r 


dr 


cos  A  (b-r) 
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A  a 
n 


lib 

-  ( — j  -  — )  /  r  sin  A  (b-r) dr 
2b  ab  a 


We  now  integrate  the  last  term  by  parts,  as  follows: 
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Using  the  eigenvalue  equation,  1.9.4,  the  expression  for  C 


n 


' 

■ 


OS 


£>  ,  nia  ft  a 


n 


J 


'■  d 


n 


:  .  oIIoS  '  j  :i  -':>X  e  !  i  ©dftipS  »-i'  WC  -  »W 

^S 


(^b 


i 

n 

aoo 

a- 

u- 

( 

0 

-  i 

n 


OS 

- -  -  0  „  .. 

»  rus  ft 


(l-d)  X  800  1  I 


:iH 


•)  (- 


•)  * 


0 


d£  dft 


t 


5  , 


n 


‘T"*“ - ,■ 

»  nis  6  -  a 


n- 


n 


»  aoo  ft  -  d 

— — .)  (- — 

A  d£  dft 


n 


<T~d)  .  nia  d 


-  ib 


n 


OS 


- ? - , - - - 

»  ale.  ft  -  d 
n 


nr.  a  -  •*>  :;oo  .6  A  -  d  A 

/jft  n  *i  **•  n  / 

( - ^ ^ - 


n 


di 

D  *xo3  noisB3'iq>zo  arid  xfr.G.I  xnoadfti/po  oul svnap 1  i  sri  •  lhibU 


4.3 


reduces  to 
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which  is  the  form  used  in  the  program. 


4 . 3  Solution  of  the  Eigenvalue  Equation 


4.2.1 


To  effect  a  numerical  solution  for  the  thermal  field 
the  roots  of  the  eigenvalue  equation,  1.9. 5, must  be  obtained. 
For  this  purpose  the  Newton-Raphson  method  proved  convenient, 
a  technique  which  will  be  found  in  any  text  on  numerical 

9 

analysis  .  The  only  difficulty  in  applying  it  to  this  par¬ 
ticular  case  was  choosing  a  starting  value  of  the  abcissa  such 
that  convergence  to  the  proper  root  was  assured. 

Examination  of  figure  1.2  shows  that  as  radius  a 
tends  to  zero  or  as  b  becomes  very  large  the  roots  approach 
the  values  nn- .  As  n  increases  the  roots  are  found  farther 
to  the  left  of  these  points,  approaching  the  values  ~^n ^ 
for  large  n.  The  problem  is  thus  to  assure  that  no  abcissa 
is  ever  specified  or  calculated  such  that  the  tangent  to  the 
curve  at  that  point  crosses  the  axis  outside  of  the  proper 
interval.  Perhaps  this  is  best  illustrated  by  an  example. 

The  method  results  in  the  algorithm 
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where  the  subscript  indicates  the  number  of  the  iteration. 
For  this  case 


ax 

f (x)  =  tan  x  +  - 

b-a 


and  hence 

df  2  a 

—  =  sec  (x)  +  - 

dx  b-a 

The  algorithm  for  this  problem  thus  becomes 
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i+1 


tan  (x . ) 
 1 

2  /  \ 
sec  (x . ) 
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ax . 
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b-a 


+ 


a 

b-a 


th 

Suppose  we  are  seeking  the  15  eigenvalue,  a,-^, 

with  a  =  50  and  b  =  300,  Further  assume  that  the  initial 

estimate,  x  ,  is  selected  as  15tt.  Then 
'  o ' 


x. 


=  15tt 


,  /-ic\  ,  50  x  15tt 

tan  ( 1 5 7T )  +  30Q  _  5Q 


sec  (15tt)  + 


50 


300  -  50 


3  TT 


1  5  TT 


1  t 


=  15tt  -  2 . 50tt 


We  know  that  a,  c  lies  between  14.5tt  and  15tt.  The 

i  5 

above  value  for  x^  is  therefore  not  acceptable,  as  further 
iteration  will  cause  convergence  not  to  but  probably  to 

a-^.  The  correction  term,  in  this  example  2.50tt,  is  much  too 
large . 
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On  the  other  hand  we  cannot  choose  the  lower  limit 

of  this  range,  xq  =  14.5tt,  because  one  of  the  denominator 

2 

terms  is  then  indeterminate  (sec  x+°°)  .  Convergence  to  the 
correct  root  will  always  occur  if  we  underestimate  it  but  stay 
above  the  particular  lower  range  limit  in  question. 

The^  computations  done  in  this  analysis  did  not  go  be¬ 
yond  twenty  terms.  Starting  values  were  chosen  one  tenth 
radian  above  the  lower  range  limit  and  proved  satisfactory 
for  all  cases. 

4 . 4  Integration 

b  sin  A  (b-r) 

The  integral  /  - - -  dr  has  been  evaluated 

a  r 

numerically  by  means  of  Simpson's  Rule.  The  interval  (a,b) 
was  divided  into  forty  parts  for  this  purpose,  applying  the 
method  twenty  times.  (Simpson's  Rule  uses  three  equally  based 
ordinates  per  application,  thus  forty  intervals  require  twenty 
applications.)  Computations  carried  out  on  this  basis  proved 
to  be  sufficiently  accurate  in  all  cases  considered. 

4 . 5  Scaling  of  Variables 

In  the  early  stages  of  the  study  the  relative  effects 
of  various  parameters  were  not  known.  For  this  reason  several 
computer  runs  were  made  varying  a  and  b  over  a  considerable 
range,  thereby  varying  the  "time  constant"  over  a  wide  range 
also.  In  order  to  render  the  program  useful  without  major 
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modification  for  all  values  of  these  variables,  computations 


were  done  on  a  scaled  basis.  Time  was  specified  as  a  fraction 

cy (b-a) ^ 


of  the  "fundamental"  time  constant,  viz.  x^ 


,  and 


atal 


radii  were  taken  as  the  tenth  points  of  the  interval  (a,b) . 

This  scheme  proved  satisfactory  for  the  cases  considered  except 
that  conversion  to  standard  time  units  was  then  necessary. 


B,  Sphere  in  an  Infinite  Medium 
4 . 6  Numerical  Difficulties 

Considerable  difficulty  was  encountered  in  programming 
the  solution  for  the  single  sphere  problem.  Two  integrations 
are  required,  the  ordinates  of  the  second  involving  the  re¬ 
sults  of  the  first.  For  this  reason  it  becomes  difficult  to 
obtain  a  high  degree  of  accuracy  as  reducing  the  base  width 
on  the  second  integral  necessitates  evaluating  the  first  at 
all  the  new  points .  This  involves  major  changes  in  the  com¬ 
puter  program. 

The  second  difficulty  lay  in  the  semi-infinite  ranges 
on  the  integrals.  In  the  interest  of  accuracy  it  was  neces¬ 
sary  to  determine  what  portion  of  the  interval  contributed 
significantly  to  the  total  area  and  to  neglect  the  remaining 
sections.  At  all  times  exponents  had  to  be  tested  to  assure 
that  they  fell  within  the  computer's  capacity. 

Experience  with  the  concentric  sphere  problem  revealed 
a  need  to  obtain  more  points  near  the  inner  sphere.  To  facil- 
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itate  this  the  radii  were  defined  as  a  geometric  rather  than 
arithmetic  sequence.  Examination  of  the  computer  output 
(Table  E.12)  shows  apparently  erroneous  results  at  the  larger 
radii.  This  may  be  attributed  to  the  increase  in  separation 
of  successive  radii,  as  is  verified  by  the  improvement  shown 
in  Table  E.13  in  which  the  integration  intervals  were  roughly 
halved.  The  last  column  of  these  tables  gives  the  final  tem¬ 
perature  for  this  radius  because  the  integrals  to  infinity  are 
assumed  identical  with  the  integrals  to  this  radius  in  the 
evaluation  of  the  transient  terms. 

Simpson's  Rule  was  also  used  for  this  problem.  The 
outer  limits  of  the  intervals  were  defined  as  a  geometric 
sequence  as  mentioned  above,  the  equal  spacing  required  by 
the  integration  technique  being  achieved  by  taking  the  aver¬ 
age  of  the  two  outer  limits  as  the  middle  abcissa. 

C.  General  Remarks 

4 . 7  Computing  Facility  and  Time  Requirements 

In  their  final  form  both  programs  were  written  in 
FORTRAN  IV  for  use  on  an  IBM  360  Model  30  computer.  Running 
time  for  each  program  is  approximately  ten  minutes  for  one 
set  of  data,  including  compiling. 
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4 . 8  Selection  of  Numerical  Values 

Since  the  problem  originally  arose  with  regard  to 

the  heating  of  soil,  parameters  typical  of  moist  soil  have 

been  used  in  the  computations.  Heat  capacity,  density,  and 

thermal  and  electrical  conductivities  were  taken  from  Carslaw 
8  3 

and  Jaeger  and  from  Rudenberg  .  The  electrode  potential  was 
selected  as  120  volts  purely  because  this  gave  a  convenient 
temperature  rise  of  75  Centigrade  degrees.  The  electrode 
radius  of  0.5  meters  was  chosen  arbitrarily. 

4 . 9  Units 

The  c.g.s.  system  of  units  was  used  because  typical 
values  were  found  specified  in  this  system.  The  radii  a  and 
b  are  in  meters  for  input-output  purposes  in  the  program,  but 
are  converted  by  the  machine  into  centimeters  for  use  in  the 
computations.  A  list  of  the  variables  is  given  in  the  table 
preceding  the  Introduction. 

4.10  Illustrations 

The  results  of  the  numerical  calculations  are  dis¬ 
played  on  the  following  graphs.  Supporting  data  may  be 
found  in  Appendix  E. 
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Radius  in  Per  Cent  of  Interval  (a,b) 
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Figure  4 . 3 
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Figure  4 . 5 

Electrode  Temperature  vs.  Time 

for  Single  Sphere  Problem 
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CHAPTER  5.  DISCUSSION  OF  RESULTS 
5 . 1  The  Meaning  of  "Time  Constant" 

The  purpose  of  this  dissertation  has  been  to  develop 
a  more  precise  prediction  of  the  thermal  behavior  of  buried 
electrodes.  Although  the  temperature-time  relation  has  been 
derived  for  the  cases  considered  it  is  perhaps  more  useful 
for  comparison  to  think  of  these  transient  fields  in  terms 
of  some  characteristic  parameter.  As  a  general  rule  the  res¬ 
ponse  of  a  body  to  a  step  increase  in  surrounding  temperature 
is  exponential  or  near  exponential  in  time.  A  convenient 
constant  for  this  sort  of  behavior  is  the  "time  constant"  for 
the  exponential  curve,  a  parameter  defined  as  the  length  of 
time  required  to  reach  the  final  temperature  at  the  initial 
rate  of  rise  or,  alternately,  as  the  time  taken  to  rise  by 
(1  -  1/e)  of  the  final  temperature  difference.  For  pure  ex¬ 
ponentials  these  definitions  are  entirely  equivalent. 

3 

Rudenberg  has  computed  the  time  constant  based  on 
the  first  definition  above.  In  this  thesis  time  constant  has 
been  interpreted  by  the  second.  The  two  results  will  not,  in 
general,  agree  because  the  characteristic  curve  is  not  a  pure 
exponential,  a  fact  which  becomes  obvious  when  a  semi-log 
plot  is  made.  A  pure  exponential  yields  a  straight  line. 


these  results  do  not. 
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5 . 2  Discussion  of  Results  for  Concentric  Spheres 

Figure  4.1  shows  a  typical  temperature-radius  relation¬ 
ship  with  time  as  parameter  for  the  concentric  spheres  problem. 
The  relative  behavior  of  the  actual  and  "fundamental"  time 
constants  (see  Section  4.5)  is  shown  in  Figure  4.2.  It  is 
apparent  from  both  these  curves  and  Figure  4.3  that  this  funda¬ 
mental  time  constant  is  not  an  accurate  or  even  consistent 
measure  of  the  actual  thermal  time  behavior. 

Figure  4.3  was  derived  in  part  from  4.2  by  the  follow¬ 
ing  method.  The  temperature-time  curve  for  the  inner  electrode 
surface  was  plotted  as  in  Figure  4.2.  From  this  graph  the 
scaled  time  at  which  the  temperature  increased  by  (1  -  1/e) 
of  its  final  difference  was  obtained  and  multiplied  by  the 
fundamental  time  constant  given  in  the  computer  output  table. 
The  real  time  so  calculated  is  plotted  in  Figure  4.3. 

The  point  (1,0)  on  Figure  4.3  is  actually  a  singu¬ 
larity  corresponding  to  coincident  spheres.  The  validity  of 
this  point  may  be  seen  by  considering  two  shells  of  nearly 
equal  radii  with  all  other  parameters  held  constant.  In 
this  arrangement  the  electric  field  is  high  and  hence  also 
the  heat  generation.  With  a  small  volume  of  material  to  heat 
the  result  is  a  very  rapid  temperature  rise,  or,  in  other 
words,  a  very  short  time  constant.  In  the  limit  as  the  shells 
become  coincident  the  time  constant  approaches  zero. 

Comparison  of  computer  results  reveals  an  interesting 
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feature.  As  long  as  the  final  temperature  and  the  ratio  of 
the  radii  b/a  are  kept  fixed  the  temperature  at  a  given  scaled 
radius  and  scaled  time  is  invariant.  For  example  the  results 
in  tables  E.l  and  E.10  are  based  on  the  same  final  temperature 
distribution  and  the  same  ratio  b/a.  Only  the  actual  magni¬ 
tudes  of  a  and  b  are  different.  The  temperature  at  r  =  a  + 

0.3 (b-a)  and  time  t  =  0.5  is  36.23°C  in  both  cases.  This 
would  suggest  the  possibility  of  drawing  up  a  table  on  a  per 
unit  basis  which  would  be  universally  applicable.  It  then 
would  be  necessary  to  compute  only  the  final  inner  electrode 
temperature  and  the  fundamental  time  constant  to  obtain  the 
transient  behavior  of  any  system  having  the  proper  ratio  b/a. 

As  previously  pointed  out  the  concentric  spheres  prob- 
lem  initially  was  undertaken  as  an  approach  to  the  single 
sphere  case.  Figure  4.3  clearly  indicates  that  the  actual 
time  constant  approaches  a  limiting  value  as  the  ratio  b/a 
increases.  The  limiting  value  shown  was  derived  from  the 
single  sphere  solution  and  appears  consistent  with  the  con¬ 
centric  sphere  curve.  It  also  illustrates  the  necessity  for 
choosing  a  ratio  b/a  of  the  order  of  40  to  50  if  an  accurate 
estimate  of  the  single  sphere  time  constant  is  to  be  obtained 
from  concentric  sphere  results. 

It  was  pointed  out  above  that,  for  a  given  ratio  b/a, 
there  is  a  consistent  relationship  between  temperature  and 
scaled  time.  This  means  that  the  actual  time  constant  is  a 
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specific  fraction  of  the  fundamental  time  constant,  a  fraction 
which  is  independent  of  the  magnitude  of  a  or  b.  If,  say,  the 
inner  radius  a  is  changed  and  b  adjusted  to  maintain  a  fixed 
ratio  b/a,  the  only  parameter  that  need  be  recalculated  is 
the  fundamental  time  constant.  This  parameter  is  defined  by 
the  relation 


That  is 


cy  (b-a) ^ 


2  ,b  ,  *  2 
cya  (-  -  1) 

Cl 


°tal 


Consideration  of  the  eigenvalue  equation,  1.9.5,  will  show 
that  the  roots  depend  only  on  the  ratio  b/a  and  not  on  the 
absolute  magnitude  of  the  inner  or  outer  radii.  This,  in 
conjunction  with  the  above  expression  for  t_,  leads  to  the 
conclusion  that  the  fundamental  (or  actual)  time  constant 
varies  directly  as  the  square  of  the  inner  radius  for  a  given 
ratio  b/a.  The  conclusion  presumably  could  be  generalized 
to  the  single  sphere  case  since  this  corresponds  to  a  limit¬ 
ing  value  of  the  ratio  b/a. 


5 . 3  Discussion  of  Results  for  Single  Sphere 

The  single  sphere  results  are  displayed  in  Figures 

4.4  and  4.5.  The  numerical  basis  used  is  given  in  Table  E-ll. 

Temperatures  of  the  electrode  surface  cannot  be  ob- 
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tained  directly  but  must  be  found  through  use  of  Figure  4.4. 
The  difficulty  lies  in  the  singularity  of  the  solution  at 
radius  a.  For  practical  purposes  this  problem  is  of  little 
consequence . 

The  time  constant  curve,  Figure  4.5,  is  plotted  by 
using  the  temperatures  derived  from  Figure  4.4  for  the  elec¬ 
trode  surface.  The  resulting  time  constant  is  shown  on  Fig¬ 
ure  4.3  and  appears  consistent  with  projections  based  on  the 
concentric  spheres. 

On  the  basis  of  calculations  using  typical  values  for 
soil  parameters  and  an  electrode  radius  of  50  centimeters, 
the  thermal  time  constant  is  of  the  order  of  38  days.  Using 
the  same  figures  the  time  constant  based  on  a  uniform  tempera¬ 
ture  rise  at  the  initial  rate  is  about  3.1  days.  The  latter 
method  therefore  appears  highly  pessimistic.  For  a  practical 
electrode,  however,  the  assumption  of  temperature  independent 
medium  parameters  must  be  considered,  as  moisture  movement 
induced  by  the  thermal  field  may  have  a  significant  influence 
on  these  parameters.  Depending  primarily  on  which  D.C.  pole 
the  electrode  forms  the  assumption  may  lead  to  an  optimistic 
or  pessimistic  value  for  the  actual  time  constant.  Estimates 
based  on  this  thesis  should  indicate  at  worst  an  order  of 
magnitude . 
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5 . 4  Generalization  of  Initial  Conditions 

For  the  purposes  of  these  analyses  the  initial  condi¬ 
tion  has  been  taken  as  constant  temperature  throughout  the 
medium.  It  is  theoretically  possible  to  use  other  conditions. 

Two  approaches  suggest  themselves.  The  first  is  to 
use  a  different  but  analytically  defined  function  such  that 
the  final  temperature  plus  this  function  yields  the  initial 
temperature  distribution.  The  solution  is  then  derived  by 
a  procedure  exactly  parallel  to  that  shown  herein. 

The  second  approach  could  be  of  great  practical  in¬ 
terest.  It  consists  of  defining  the  initial  condition  numer¬ 
ically  instead  of  analytically  and  performing  the  integrations 
using,  say,  tabulated  values.  This  could  be  employed  to  ad¬ 
vantage  in  cases  where  soil  parameters  change.  For  instance 
the  conditions  specified  in  this  thesis  could  be  used  to  find 
the  temperature  at  some  time  t.  Then,  taking  the  temperature 
so  derived  as  an  initial  condition  and  suitably  altering  the 
"constants",  the  temperature  at  some  later  time  could  be  de¬ 
rived.  By  choosing  fine  enough  intervals  of  time  an  excellent 
prediction  of  actual  behavior  should  be  possible.  This  tech¬ 
nique,  however,  does  not  remove  the  restrictions  imposed  by 
the  assumption  of  homogeneity. 

The  same  analysis  can  be  used  to  predict  cooling  of 
the  electrode  system.  It  is  only  necessary  to  set  the  final 
temperature  (steady  state)  to  a  constant  and  reverse  the  sign 
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of  the  transient  portion  of  the  solution.  This  could  be  used 
to  predict  the  necessary  inter-use  intervals  if  intermittent 
operation  was  necessary. 

These  possibilities  apply  only  if  the  same  boundary 
conditions  are  to  hold  for  the  auxiliary  problems.  If  this 
is  not  the  case  almost  the  whole  analysis  must  be  repeated. 

In  that  event  this  should  serve  as  a  guide  since  the  essential 
ideas  are  the  same. 

5 . 5  Finite  Cylinder 

The  solution  for  a  finite  cylinder  was  attempted  but 
the  mathematical  difficulties  prevented  any  progress.  For 
practical  purposes,  however,  the  solutions  derived  in  this 
thesis  may  be  of  some  use.  A  little  thought  will  show  that, 
for  regions  very  close  to  a  long  but  finite  rod,  the  thermal 
field  must  behave  like  that  of  the  infinite  cylinder.  For 
regions  well  removed  from  the  rod  the  situation  is  similar  to 
that  of  the  sphere  in  an  infinite  medium. 

5 . 6  Laplace  and  Fourier  Transform  Techniques 

Two  other  approaches  to  the  problems  solved  herein 
are  open,  viz.  Laplace  and  Fourier  Transforms.  However,  all 
attempts  at  using  them  led  to  very  difficult  complex  plane 
integrals  and  hence  they  were  abandoned  in  favor  of  the  tech¬ 
nique  of  reducing  the  problem  to  a  homogeneous  case  by  trans¬ 


formations  . 
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Miller‘S  shows  an  approach  for  the  extension  of  a 
bounded  region  problem  to  an  infinite  medium  situation  by 
working  from  a  standard  trigonometric  Fourier  series  solution 
into  the  Fourier  integral.  This  approach  was  attempted  (by 
Professor  Berg)  but  difficulties  arose  in  applying  the  boun¬ 
dary  condition  at  the  inner  surface.  The  attempt  was  abanr 
doned  when  success  was  obtained  with  the  method  shown  herein. 

5 . 7  General 

Although  specific  reference  is  made  to  soil  electrodes 
the  results  apply  to  any  medium  which  meets  the  assumed  con¬ 
ditions.  Application  to  A.C.  is  undoubtedly  possible  pro¬ 
viding  the  frequency  is  not  too  high.  If  generalization  to 
this  area  is  contemplated  it  would  be  advisable  to  review  the 
entire  analysis  in  that  light. 
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CHAPTER  6.  SUMMARY  AND  SUGGESTIONS 
FOR  FURTHER  DEVELOPMENT 

6 . 1  Concentric  Spheres 

The  analysis  of  Chapter  1  resulted  in  the  following 
equations  for  concentric  spheres: 


2 

-A  kt 
n 

oo  e 

T(r,t)  =  l  - 

n=l  r 


C  sin  A  (b-r) 
n  n 


where 


G  G 


G  G 


- T  -  +  T  ,  +  - TT  ~  - 
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2r  ar  2b  ab 
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b-a  sin  A  (b-a)  a 
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/  r  sin  A  (b-r) 
j  n 


2r 


G  G  G 

- 2  +  — 

ar  2b  ab 


dr , 


o  U  ab  „ 

G  =  S.  (-2_ )2 


°t  b-a 


th 


and  An  is  the  n  root  of  the  equation 


A  a  cos  A  (b-a)  +  sin  A  (b-a)  =  0 

n  n  n 

The  remaining  symbols  are  summarized  (page  vi)  in  a 
table  immediately  preceding  Chapter  1.  A  computer  program 
of  this  solution  may  be  found  in  Appendix  D. 
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6 . 2  Sphere  in  an  Infinite  Medium 


The  following  equations  were  developed  in  Chapter  2.: 
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Appendix  D  contains  a  computer  program  of  this  solution. 


6 . 3  Concentric  Infinitely  Long  Cylinders 


Chapter  3  contains  the  development  of  the  following 

equations  describing  the  transient  thermal  behavior  of  the 
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region  between  two  concentric  infinitely  long  cylinders. 
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and  A^  satisfies  the  equation 
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No  computer  program  or  numerical  solution  is  given. 


6 . 4  Suggestions  for  Further  Development 

Perhaps  the  most  important  problem  resulting  from 
this  analysis  is  that  of  experimental  verification.  The  con¬ 
centric  sphere  problem  should  lend  itself  quite  nicely  to 
this  task. 

Assuming  a  positive  result  to  such  an  endeavor  it 
would  be  of  considerable  practical  interest  to  consider  the 
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validity  of  the  equivalent  sphere  relations  given  by  Rudenberg 
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with  regard  to  transient  thermal  behavior  or  to  develop  more 

2 

accurate  expressions  for  design  purposes.  Taylor's  sug¬ 
gestion  of  using  "specific  loading"  may  be  useful  in  this  res¬ 
pect. 

The  possibility  of  summarizing  the  results  of  this 
thesis  in  tabular  form  as  discussed  in  Section  5,2  may  be 
tempting  from  the  point  of  view  of  operation  of  electrodes. 
Closely  tied  with  this  aspect  is  the  numerical  solution  for 
different  initial  conditions  as  discussed  in  Section  5.4. 
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APPENDIX  A.  DEVELOPMENT  OF  EQUATIONS 


A , 1  Electric  Equations 


Start  with  one  of  Maxwell's  equations 
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and  Ohm's  law 
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From  A. 1.1  and  A. 1.2  the  result 
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follows . 

Using  A. 1.3  and  A. 1.4  and  expanding  this  by  the  chain 
rule  of  differential  calculus  gives 
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or , rearranging  the  terms 
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For  steady  state  fields 
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Thus  A. 1.5  becomes 
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which  is  equivalent  to,  from  A. 1.4 
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If  we  assume  a  is  invariant  A. 1.6  reduces  to 
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or,  in  terms  of  potential 
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A. 2  Equivalence  of  Specifying  U  or  I  in  the 
Heat  Generation  Term 


Ohm's  law  states  that 
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that  is  the  total  electrode  current  must  flow  through  any 
closed  surface  surrounding  the  electrode.  At  a  fixed  radius, 
considering  the  radial  symmetry  of  the  cases  in  question,  the 
quantity  agE  is  constant  in  magnitude  and  always  parallel  to 
the  surface  element  dS .  Using  this  fact  and  equation  A. 2. 2 
in  A. 2.1,  it  follows  that 
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Consider  first  the  concentric  spheres 
The  area  of  a  sphere  is  given  by 
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Now,  from  equation  1.4.2 
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The  current  flowing  from  the  inner  electrode  for  the  concentric 
sphere  case  is  thus  given  by 
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The  constant  G  may  be  expressed  in  terms  of  current 
by  the  following  procedure.  G  was  defined  in  Section  1.4  by 
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Substituting  from  equation  A. 2. 4  this  may  be  written  as 
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that  is  in  terms  of  current  rather  than  potential. 

Consider  now  the  single  sphere.  Again  the  area  is 

given  by 
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The  electrode  current  is  therefore 
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G  for  this  case  was  defined  in  Section  2.3  by 
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which  becomes,  using  the  above  relation 


/  2 
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The  equivalence  of  specifying  Uq  or  I  is  thus  estab¬ 
lished.  For  practical  purposes  the  current  I  is  probably 
more  useful. 
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Differentiating  with  respect  to  t  yields 
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which,  when  compared  with  A. 3.1,  will  be  found  identical. 
The  desired  proof  is  thus  effected. 


A . 4  Derivation  of  the  Diffusion  Equation 


Fourier's  equation  states  that  the  rate  of  heat  trans¬ 
fer  is  proportional  to  the  gradient  of  temperature.  Mathe¬ 
matically 
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The  rate  of  heat  generated  per  unit  volume  by  the  passage 
of  an  electric  current  is  given  by 
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and  the  rate  of  heat  stored  per  unit  volume  by 
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Consider  a  rectangular  parallepiped  as  shown  in  Figure  A.l. 
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Figure  A.l 

Differential  Element  for  the  Development 

of  the  Diffusion  Equation 
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The  energy  balance  for  this  elemental  volume  may  be  written 
as  follows: 
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We  evaluate  the  partial  space  derivatives  by  means  of  Taylor' 
Theorem,  neglecting  terms  of  order  greater  than  two.  This 
results  in  the  following  six  equations. 
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If  we  recognize  the  space  dependence  as  the  operator  V  this 
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The  term  on  the  right  hand  side  is  normally  given  as 
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is  called  the  thermal  diffusivity  of  the  medium.  The  thermal 
field  equation  is  thus 
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APPENDIX  B.  ORTHOGONALITY  OF  SIN  A  (b-r)  OVER  (a,b) 


We  wish  to  show  that 
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We  begin  by  integrating  by  parts.  Let 
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and 
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Comparing  the  integral  on  the  right  with  the  original,  it 
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Substituting  this  in  the  above  relation  gives 
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From  the  eigenvalue  equation  we  have 
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which,  when  employed  in  the  above  relation  gives 
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This  is  the  required  orthogonality  constant.  We  have  thus 
established  the  desired  relations,  B.1.1. 
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APPENDIX  C.  ORTHOGONALITY  OF  W  OVER  (a,b) 
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We  wish  to  evaluate  the  orthogonality  constant  and  prove 
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Following  the  process  in  Carslaw  and  Jaeger  ,  we  know  that 
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Integrating  both  sides 
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and  we  know  from  the  eigenvalue  equation,  3.9.8  that 
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C  HEAT  FLOW  IN  CONCENTRIC  SPHERES.  A.  HERRON  FOR  G.  BERG 

DIMENSION  ALPHA (20) ,  TRT( 11),TRF(11),Y(21,20),C(20) 

53  READ(1,1  )SIGE,SIGT,U0,R1,R0,Z,GAM,TAMB 

1  FORMATf E10.2,7F10.6) 

WRITE(3,67) 

67  FORMAT ( 1H 1 ) 

WRITE(3,2) 

2  FORMAT ( 1H9 , 19H  SIGMA  EL  SIGMA  TH , 5X , 2HU0 , 8X , 2HR1 , 8X , 2HR0 , 9X , 1HC , 

16X,5HGAMMA,6X,5HT  AMB,5X,5HR1/R0,7X,1HI) 

2,  /  R  0 

F  1=4  .*3 . 14159  3*U0*R1*R0  /  (Rl-RO  )‘*SIGE 
WRITE ( 3,50)SIGE,SIGT,U0,R1 , RO , Z , GAM , T AMB , RAT , FI 
50  FORMAT ( 1H0,E10.3,8F10.4,E10.3) 

R1=R1*100. 

R0=R0*100. 

WRITE ( 3,3 ) 

3  FORMAT ( 1H0 ,9HALPHA  N'S) 

B  =  RO/ (Rl-RO ) 

DO  20  1=1,20 
FI= 1*2-1 

ALPHA ( I )=F I *1.57079 63+0.1 
21  TEST=ALPHA ( I ) 

D=1 • /COS ( ALPHA ( I ) ) 

TAN  =  S IN ( ALPHA ( I  )  ) /COS ( ALPHA ( I ) ) 

ALPHA ( I ) =AL PHA ( I  )-(  ( TAN+B*AL PHA ( I ) )/(B  +  D*D)  ) 

X= ABS ( 1 .-TEST /ALPHA ( I ) ) 

I F ( X— 1 • E— 03 ) 20 , 20 , 2 1 
20  CONTINUE 

WRITE ( 3 ,4) ALPHA 

4  FORMAT ( 1H0, 10F10. 3 ) 

TAU=Z*GAM* (Rl-RO ) * (Rl-RO ) / ( SI GT* ALPHA ( 1 ) *  A  L  P  H  A ( 1 )  )/ 86400 . 

WRITEt 3,5  )TAU 

5  FORMAT ( 12H0TAU  ( DA YS ) = , E 14 . 7 ) 

WRITE( 3,6) 

6  FORMAT ( 1 HO , 6H  T I  ME , 5X , 1H0 , 6X , 4H1 / 1 0 , 4X , 4H2 / 1 0 , 4X , 4H3 / 1 0 , 4X , 
14H4/10,4X,4H5/10,4X,4H6/10,4X,4H7/10,4X,4H8/10 , 4X , 4H9 / 1 0 , 5X , 1H1 ) 

V4=S I GE*U0*U0*R1*R 1*R0*R0/ ( SIGT* (Rl-RO)* ( Rl-RO ) *4. 1852 ) 

CONST  =  T AMB  +  V4* (  ( 1 . / ( 2 .*R 1*R 1 ) ) - ( 1 . / ( RO*R 1 )  )  ) 

DO  22  J=l,ll 
GJ= J-l 
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FJ=RO+GJ* ( Rl-RO ) /10. 

22  TRF( J)=V4*(  ( 1  .  / ( R  0  *  F  J )  )-( 1 . / ( 2 .*F J*F J )  )  )+CONST 
DO  25  K=  1  ,  13 
IF(K-5)61,61,62 

61  T IME=K-1 

T I  ME  =  T  I  ME/40 • 

GO  TO  63 

62  T I ME=K-3 

T IME=T I ME*0 .05 

63  DO  25  L= 1 , 20 

WATCH=ALPHA(L )*ALPHA( L ) /ALPHA ( 1 ) /ALPHA ( 1 ) ^  T I M  E 
IF (WATCH-99. ) 5 1,52 ,52 
52  Y(K,L)=0. 

GO  TO  25 

51  Y(K,L )=EXP(-WATCH) 

25  CONTINUE 

H= (Rl-RO)/ 40 . 

DO  31  N= 1 , 20 
SUM=0. 

R  =  RO 

DO  30  M= 1 ,39,2 
F  M  =  M  - 1 

VNT1=SIN( ( 1.-FM/40. )*ALPHA(N) )/R 
R  =  R  +  H 

VNT2=S IN ( ( 1 . - ( FM+1 • )/40. )*ALPHA(N) )/R*4. 

R  =  R+H 

VNT3=SIN( ( 1 .- ( FM+2 • )/40. )*ALPHA(N) )/R 

30  SUM=SUM+ ( VNT1+VNT2+VNT3 )*H/3. 

C0EFC=V4*2./ (R1-R0*SIN( ALPHA(N) ) *S IN ( ALPHA ( N ) ) ) 
TERM2= ( COS  (ALPHA(N) )-l. ) / ( R0*AL PH A ( N ) )*( Rl-RO) 
TERM3=( Rl-RO )*{ 1 . /RO- 1 . / ( 2 . *R 1 ) ) /ALPHA ( N ) 

31  C  (  N  )  =COEFC*  (  SUM/2. +TERM2+TERM3  ) 

TRT( 11 ) =TRF ( 11 ) 

DO  32  I T  =  1 , 13 
IF (  IT-5 )64,64,65 

64  FIT  = I T— 1 
FIT=FIT/40. 

GO  TO  66 

65  F I T= I T-3 
FIT=FIT*0.05 
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DO  33  IR=1,10 
F I R= I R— 1 

ss=o. 

DO  34  IN=1,20 

S  =  Y(  IT, IN)*C ( IN)*SIN  ( ALPHA ( IN) 4(1. -FIR/10.)  ) 
SST=SS 
SS=S+SS 

ST=ABS  ( l.-SST/SS ) 

IF ( ST-l.E-03)37,37,34 
34  CONTINUE 

37  SS=SS/(R0+FIR*(R1-R0)/10. ) 

33  TRT ( IR)=SS+TRF( IR) 

32  WRITE(3,7)FIT,TRT 

7  FORMAT { 1H0»F7.4,11F8.2) 

8  FORM AT ( 8H0  INF  IN , 1 1F8 .2// ) 

60  FORMAT { 13) 

WR I T E ( 3  »  8 ) TRF 
RE AD ( 1 ,60  )NC 
IF  ( NC ) 54, 53  »  53 
54  CALL  EXIT 

END 
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C  HEAT  FLOW  FOR  SPHERE  IN  INFINITE  MEDIUM  A.  HERRON  FOR  G.  BERG 

DIMENSION  R(22),TEMP(9  ),FT(9  ) , V ( 1 1 ) , PH  I ( 22 ) , SUM ( 22  ) 

12  RE AD ( 1 , 1 )SIGE,SIGT,UO,RO,HC,GAM,TAMB,TAU 

1  FORMAT (E10.2,6F10.6,E10.3) 

WR  ITE ( 3  ,  360  ) 

360  FORMAT ( 1H 1  ) 

WRITE (3, 2) 

2  FORMAT ( 20H9  SIGMA  EL  SIGMA  TH , 5X , 2HU0 , 8X , 2HR0 , 8X , 2HHC , 

16X, 5HGAMMA,6X,5HT  AMB,7X, 1HI/  ) 

F I =4. *3, 141593*SIGE*U0*R0 

WRITE ( 3,3 )SIGE,SIGT,UO,RO,HC,GAM,TAMB, FI 

3  FORMAT ( 1H  , E 10 . 4 , 7F 10 . 4/ / ) 

DO  4  1=2,18,2 

11=1-2 

4  R (  I ) =R0  *( 1. +  .1*1. 414214**1 1 ) 

WRITE ( 3,5  )  (R ( I  )  ,  1=2,18,2) 

5  FORMAT ( 1 1 H  TIME/RAD  ,  9E10.4//) 

R ( 1 )  =  ( RO  +  R ( 2  )  )  * 0 . 5 

DO  225  1=3,17,2 

225  R( I )  =  ( R( 1  +  1 )+R( 1-1  )  )*0.5 
R0=R0*100 . 

DO  203  1=1,18 
203  R(I)=R(I)*100. 

CAY3=S IGT/HC/GAM 

CAY4=S I GE*UO*UO*RO*RO/S IGT / 4. 1852 

DO  6  J= 1 , 9 

JJ=2*J 

6  FT(J)=-CAY4/2./R ( J J ) +TAMB*R ( J J ) +CA Y4/R0 
DO  7  LT= 1 ,9 

I F ( LT-1 )356,356,357 

356  T I ME=0 . 1*TAU 
GO  TO  358 

357  TIME=LT-1 
TIME=TIME*TAU/4. 

358  COEF 1=C A Y4*R0*0. 5/1 .772454 
COE  F2  =  CAY4=:;0. 5/ 1.772454 
C0EF3=-CAY4/R0/1 .772454 
DEN2=0.25/CAY3/TIME 

DEN  =  SQRT (DEN 2  ) 

DO  316  LR= 1 ,18 
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351 


G  =  0. 

ARE A=0 . 

INT=-1 

H=( (R(LR)-RO)*(R(LR)-RO) )*DEN2 
IF ( H-90. ) 321 , 308 , 308 

321  H= ( SORT (90. ) / DEN+RO-R ( LR ) )/40. 

309  FUNC=COEF 1 / ( ( G+RO ) * ( G+RO ) ) +C0EF2/ ( G+RO ) +C0EF3 
POWERS ( (R(LR)-RO+G)*(R(LR)-RO+G) )*DEN2 

IF ( POWER-90 • ) 300 , 301 , 301 

300  ETERM=EXP ( -POWER ) 

GO  TO  302 

301  ETERM=0. 

302  IF ( INT ) 303 ,304, 305 

303  FB=FUNC*ETERM 
G  =  G+H 

I  NT  =  0 
GO  TO  309 

304  FM=FUNC*ETERM 
G  =  G  +  H 

I NT=+ 1 
GO  TO  309 

305  FE=FUNC*ETERR 

ARE A= ARE A- ( FB+4. *FM+FE ) *H/3 . 

307  FB  =  F  E 
I  NT  =  0 
G  =  G+H 

P 0 W E R  =  (  (R(LR)-RO+G)*(R(LR)-RO+G)  )*DEN2 
I  F ( POWER -9  0. ) 309 ,308,3  08 

308  G=R(LR)-R0-3.*SQRT( 11. ) /DEN 
IF ( G) 322 ,323,323 

322  H=(R(LR)-R0+3.*SQRT( 11. )/DEN)/40. 

G=0 . 

GO  TO  324 

323  H=6.0*SQRT( 11. J/DEN/40. 

324  I NT=- 1 

310  FUNC  =  CO E F 1 / (  ( G  +  RO ) *( G+RO  )  ) +C0EF2/ ( G  +  RO ) +C0EF3 
POWER= ( ( R ( LR )-RO-G )* ( R ( LR )-RO-G ) )*DEN2 

I F ( POWER-9  0 . ) 353 ,354,354 
354  ETERM=0 . 

GO  TO  355 
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353  ETERM  =  EXP (-POWER  ) 

355  I  F (  INT)311,312,313 

311  FB=FUNC*ETERM 
G  =  G  +  H 

I  NT  =  0 
GO  TO  310 

312  FM=FUNC4ETERM 
G  =  G  +  H 

I  NT=  1 
GO  TO  310 

313  FE=FUNC*ETERM 

ARE A= ARE A+ ( FB+4.4FM+FE ) 4H/3. 

315  FB=FE 
I  NT  =  0 
G  =  G  +  H 

POWER=(  (R(LR)-RO  — GHMR(LR)-RO-G)  J4DEN2 
IF ( POWER-90. )310,316,316 

316  PHI ( LR)=AREA*DEN 

219  SUM ( 2 ) = ( 4. 4PH I ( 1 ) 4EXP ( -R ( 1 ) /RO ) +PH I ( 2 ) 4EXP ( -R ( 2 ) / RO ) ) 4 ( R ( 2 ) -R ( 1 ) ) / 

13. 

DO  101  M=2,16,2 
ARG1=PHI ( M ) 4EXP ( -R ( M ) /RO ) 

ARG2=PHI (M  +  l )*EXP (-R(M+1 ) /RO ) 

IF(R(M+2)/R0-99. )112,112,210 

210  ARG3  =  0 . 

GO  TO  101 

112  ARG3  =  PH  I (M  +  2 )4EXP ( -R (M  +  2 ) /RO  ) 

101  SUM (M+2 )=SUM (M )+(ARG 1+4.4 ARG2+ARG3 ) / 3 . 4 ( R ( M+l ) -R ( M ) ) 

221  DO  205  LR=1,9 
LRR=24LR 

IF(R(LRR)/RO— 99. )211,211,212 
212  V ( L  R ) =0 . 

GO  TO  205 

211  V ( LR ) = ( SUM ( 1 8 ) -SUM ( LRR ) ) /R04EXP ( R ( LRR ) /RO ) 

205  TEMP(LR)=(FT(LR)+V(LR) ) /R ( LRR ) 

TIME  =  T  I  ME/86400. 

7  WRITE(3,9)TIME,TEMP 

DO  202  M= 1 » 9 
M  M = 2  4  M 

202  FT(M)=FT(M)/R(MM) 
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FORMAT ( 1H0,E10.4,11F10.2) 

FORMAT ( 11 HO  INFINITY  ,  1  IF  10. 2//  ) 
WRITE ( 3, 10 ) FT 
RE AD ( 1,359 )NC 
FORMAT (  13  ) 

I F ( NC ) 13 , 12 , 1 2 
CALL  EXIT 
END 
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C  IMPROVED  ACCURACY 

C  HEAT  FLOW  FOR  SPHERE  IN  INFINITE  MEDIUM  A.  HERRON  FOR  G.  BERG 

DIMENSION  R(36),TEMP(9  ),FT(9  )  , V ( 1 1 ) , PH  I ( 36 ) , SUM ( 36 ) 

12  READ! 1,1 )SIGE,SIGT,UO,RO,HC,GAM,TAMB,TAU 

1  FORMAT (E10.2,6F10.6,E10.3) 

WR I TE ( 3 , 360 ) 

360  FORMAT ( 1H1 ) 

WR  ITE ( 3 , 2 ) 

2  FORMAT ( 20H9  SIGMA  EL  SIGMA  TH , 5 X , 2HU0 , 8X , 2HR0 , 8X , 2HHC , 
16X,5HGAMMA,6X,5HT  AMB,7X,1HI/) 

FI =4. *3. 141593*SIGE*U0*R0 

WRITE ( 3,3)SIGE,SIGT,U0,R0,HC,GAM,TAMB,FI 

3  FORMAT ( 1H  ,  E 1 0 . 4 , 7F 1 0 . 4/ / ) 

DO  4  1=2,36,2 

11  =  1-4 

4  R (  I  )  =R0  *( l.  +  . 1*1. 189207**1  I ) 

WR I TE ( 3 , 5 ) ( R ( I ) ,1=4,36,4) 

5  FORMAT ( 1 1H  TIME/RAD  ,  9E10.4//) 

R( 1 )=(R0+R(2) 

DO  225  1=3,35,2 

225  R(I)=(R(I+1)+R(I-1))*0.5 

R0=R0*100. 

DO  203  1=1,36 
203  R ( I ) =R ( I ) *100 • 

CAY3=SIGT/HC/GAM 

CAY4=SIGE*U0*U0*R0*R0/SIGT/4. 1852 
DO  6  J= 1 , 9 
J J=4* J 

6  FT( J)=-CAY4/2./R( JJ)+TAMB*R( JJ)+CAY4/R0 
DO  7  LT= 1 ,9 

IF(LT-1)356,356,357 

356  TIME=0.1*TAU 
GO  TO  358 

357  T I ME=LT-1 
TIME=TIME*TAU/4. 

358  C0EF1=CAY4*R0*0. 5/1 .772454 
C0EF2=CAY4*0. 5/ 1.772454 
C0EF3=-CAY4/R0/ 1.772454 
DEN2=0.25/CAY3/TIME 
DEN=SQRT ( DEN2 ) 

DO  316  LR=1 »  36 
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351 


G=0. 

ARE A  =  0 . 

I N  T  =  - 1 

H=(  (R(LR)-RO)*(R(LR)-RO)  )*DEN2 
IF  (IH-90.  >321,308,308 

321  H= ( SORT ( 9  0 . )/DEM  +  RO-R( LR)  )/40. 

309  F  U  N  C  =  C  0  E  F 1 / (  ( G+RO ) *< G+RO )  ) +COE F2/ ( G+RO ) +C0EF3 
POWER= ( (R(LR)-RO+G)*(R(LR)-RO+G) )*DEN2 

IF ( POWER-90 •  )  300 ? 301 , 301 

300  ETERM=EXP( -POWER) 

GO  TO  302 

301  E  T  E  R  M  =  0 . 

302  IF(  INT)303,304,305 

303  FB=FUNC*ETERM 
G=G  +  H 

I  N  T  =  0 
GO  TO  309 

304  FM=FUNC*ETERM 
G=G  +  H 

I NT=+ 1 
GO  TO  309 

305  FE=FUNC*ETERM 

ARE A= ARE A- (FB+4.*FM+FE)*H/3. 

307  FB  =  F  E 
I  NT  =  0 
G-G  +  H 

POWER=  (  (R(LR)-RO+G)*(R(LR)-RO+G)  )*DEN2 
IF ( POWER— 90. >309,308,308 

308  G=R(LR)-R0-3.*SQRT( 11. ) /DEM 
IF(G)322,323,323 

322  H=(R( LR)-R0+3.*SQRT ( ll.)/DEN)/40. 

G  =  0  • 

GO  TO  324 

323  H=6.0*SQRT( 11. )/DEN/40. 

324  INT=-1 

310  FUNC  =  C0EF1 / (  ( G  +  RO ) * ( G+RO  )  ) +C0EF2/ ( G  +  RO ) +C0EF3 
POWER=  (  (R(LR)-RO-G)*<R(LR)-RO-G)  )*DEN2 

IF ( POWER-90. >353,354,354 
354  ETERM=0. 

GO  TO  355 
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353  ETERM=EXP( -POWER) 

355  I F ( INT)311, 312,313 

311  FB=FUNC*ETERM 
G  =  G+H 

I  NT  =  0 
GO  TO  310 

312  FM=FUNC*ETERM 
G  =  G  +  H 

I  NT=  1 
GO  TO  310 

313  FE=FUNC*ETERM 
AREA=AREA+(FB+4.*FM+FE )*H/3. 

315  FB=FE 
I  NT  =  0 
G=G+H 

POWER= ( (R(LR)-RO-G)*(R(LR)-RO-G) )*DEN2 
I  F ( POWER -9  0. )310,316,316 

316  PHI ( LR )=AREA*DEN 

219  SUM(2)=(4.*PHI (1)*EXP(-R( 1)/R0)+PHI (2 )*EXP (-R ( 2 ) /RO ) ) * ( R ( 2 ) -R ( 1 )  )  / 

13. 

DO  101  M=2,34,2 
ARG1=PHI (M)*EXP(-R(M)/RO) 

ARG2=PH I (M+1)*EXP(-R(M+1 )/R0) 

IF ( R ( M+2 ) /RO-99. )112,112,210 

210  ARG3=0 • 

GO  TO  101 

112  ARG3  =  PH  I  (M  +  2  )  *  E  X  P  (  -R  (  M  +  2  )  /RO  ) 

101  SUM ( M  +  2  )  =SUM ( M )  +  ( ARG 1+4 . * ARG2  +  ARG3 ) / 3 . * ( R ( M+l )-R(M) ) 

221  DO  205  LR=1,9 
LRR=4*LR 

IF (R(LRR) /RO-99. )211,211,212 
212  V(LR)=0. 

GO  TO  205 

211  V(LR)=(SUM(36)-SUM(LRR) ) /RO*EXP ( R ( LRR ) /RO ) 

205  TEMP(LR)=(FT(LR)+V(LR) )/R(LRR) 

T I ME  =  T I  ME / 8 6400 • 

7  WRITE(3,9 )TIME,TEMP 

DO  202  M=l,9 
MM=4*M 

202  FT (M)=FT { M ) /R ( MM ) 
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9  FORMAT ( 1H0,E10.4, 11F10.2 ) 

10  FORMAT  (  11 F10  INFINITY  ,  11F10.2//) 
WRITE ( 3, 10)FT 

READ ( 1 ,359 )NC 
359  FORMAT (13) 

IF(NC)  13,12,12 
13  CALL  EXIT 

END 
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APPENDIX  E 


NUMERICAL  TABLES 

Input-Output  Variables 

Symbol 

Parameter 

Units 

SIGMA  EL 

electrical  conductivity 

1/ (ohm-cm 

SIGMA  TH 

thermal  conductivity 

cal/  ( sec-' 

UO 

electric  potential 

volts 

Rl 

outer  shell  radius 

m . 

RO 

inner  shell  radius 

m. 

C ,  HC 

heat  capacity 

cal/gm 

GAMMA 

mass  density 

gm/cm^ 

Tamb 

ambient  temperature 

°c 

I 

electrode  current 

amperes 

Tau 

fundamental  time  constant 

days 
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Table  E.l  Concentric  Spheres 
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